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We formulate low Mach number fluctuating hydrodynamic equations appropriate for mod- 
eling diffusive mixing in isothermal mixtures of fluids with different density and transport 
coefhcients. These equations eliminate the fast isentropic fluctuations in pressure associ- 
ated with the propagation of sound waves by replacing the equation of state with a local 
thermodynamic constraint. We demonstrate that the low Mach number model preserves the 
spatio-temporal spectrum of the slower diffusive fluctuations. We develop a strictly conserva- 
tive finite- volume spatial discretization of the low Mach number fluctuating equations in both 
two and three dimensions. We construct several explicit Runge-Kutta temporal integrators 
that strictly maintain the equation of state constraint. The resulting spatio-temporal dis- 
cretization is second-order accurate deterministically and maintains fluctuation-dissipation 
balance in the linearized stochastic equations. We apply our algorithms to model the devel- 
opment of giant concentration fluctuations in the presence of concentration gradients, and 
investigate the validity of common simpliflcations neglecting the spatial non-homogeneity of 
density and transport properties. We perform simulations of diffusive mixing of two fluids of 
different densities in two dimensions and compare the results of low Mach number continuum 
simulations to hard-disk molecular dynamics simulations. Excellent agreement is observed 
between the particle and continuum simulations of giant fluctuations during time-dependent 
diffusive mixing. 
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I. INTRODUCTION 

Stochastic fluctuations are intrinsic to fluid dynamics because fluids are composed of molecules 
whose positions and velocities are random at thermodynamic scales. Because they span the whole 
range of scales from the microscopic to the macroscopic [U |2] , fluctuations need to be consistently 
included in all levels of description. Stochastic effects are important for flows in new microfluidic, 
nanofluidic and microelectromechanical devices P] ; novel materials such as nanofluids [1] ; biological 
systems such as lipid membranes [5J, Brownian molecular motors [6j, nanopores [7]; as well as 
processes where the effect of fluctuations is amplified by strong non-equilibrium effects, such as 
ultra clean combustion, capillary dynamics [HI [9], hydrodynamic instabilities |10H12j . and others. 

One can capture thermal fluctuations using direct particle level calculations. But even coarse- 
grained particle methods [Tl[13l[Tl] are computationally expensive because the dynamics of individ- 
ual particles is significantly faster than hydrodynamic time scales. Alternatively, thermal fluctua- 
tions can be included in the Navier-Stokes equations through stochastic forcing terms, as proposed 
by Landau and Lifshitz [15] and later extended to fluid mixtures [16] . The basic idea of fluctuating 
hydrodynamics is to add a stochastic flux corresponding to each dissipative (irreversible, diffusive) 
flux |17j . This ensures that the microscopic conservation laws and thermodynamic principles are 
obeyed while also maintaining fluctuation-dissipation balance. Specifically, the equilibrium thermal 
fluctuations have the Gibbs-Boltzmann distribution dictated by statistical mechanics. Fluctuating 
hydrodynamics has been demonstrated to be a very useful tool in understanding complex fluid flows 
far from equilibrium [16]. Theoretical calculations are often very tedious and only feasible after 
ignoring nonlinearities, inhomogeneities in density, temperature, and transport properties, surface 
dynamics, gravity, unsteady flow patterns, and other important effects. In the past decade fluc- 
tuating hydrodynamics has been applied to study rather nontrivial practical problems [9l [T8H20] ; 
however, the numerical methods used are far from the state-of-the-art in deterministic solvers. 

Previous computational studies of the effect of thermal fluctuations in fluid mixtures [3 HH 
[2T] have been based on the compressible fluid equations and thus require small time steps to 
resolve fast sound waves (pressure fluctuations). Recently, some of us developed finite- volume 
methods for the incompressible equations of fluctuating hydrodynamics [22J, which eliminate the 
stiffness arising from the separation of scales between the acoustic and vortical modes [23^ I24j . 
For inhomogeneous fluids with non-constant density, diffusive mass and heat fluxes create local 
expansion and contraction of the fluid and the incompressibility constraint should be replaced by a 
"quasi- incompressibility" constraint [24:\ [25] . The resulting low-Mach number equations have been 
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used for some time to model deterministic flows with thermo-chemical effects [241 126j , and several 
conservative finite- volume techniques have been developed for solving equations of this type [271 - 
131] ■ To our knowledge, thermal fluctuations have not yet been incorporated in low Mach number 
models. 

In this work we extend the staggered-grid, flnite-volume approach developed in Ref. |22j to 
isothermal mixtures of fluids with unequal densities. The imposition of the quasi-incompressibility 
constraint poses several nontrivial mathematical and computational challenges. At the mathemat- 
ical level, the traditional low Mach number asymptotic expansions \23\ [23] assume spatio-temporal 
smoothness of the flow and thus do not directly apply in the stochastic context. At the com- 
putational level, enforcing the quasi-incompressibility or equation of state (EOS) constraint in a 
conservative and stable manner requires specialized spatio-temporal discretizations. By careful se- 
lection of the analytical form of the EOS constraint and the spatial discretization of the advective 
fluxes we are able to maintain strict local conservation and enforce the EOS to within numerical 
tolerances. In the present work, we employ an explicit projection-method temporal discretizations 
because of the substantial complexity of designing and implementing semi-implicit discretizations 
of the momentum equation for spatially-inhomogeneous fluids. 

Thermal fluctuations exhibit unusual features in systems out of thermodynamic equilibrium. 
Notably, external gradients can lead to enhancement of thermal fluctuations and to long-range 
correlations between fluctuations |16l [32^135] . Sharp concentration gradients present during diffu- 
sive mixing lead to the development of macroscopic or giant fluctuations |36H38] in concentration, 
which have been observed using light scattering and shadowgraphy techniques (21 [391 HO] ■ These 
experimental studies have found good but imperfect agreement between the predictions of a sim- 
plifled fluctuating hydrodynamic theory and experiments. Computer simulations are, in principle, 
an ideal tool for studying such complex time-dependent processes in the presence of nontrivial 
boundary conditions without making the sort of approximations necessary for analytical calcula- 
tions, such as assuming spatially-constant density and transport coefficients and spatially-uniform 
gradients. On the other hand, the multiscale (more precisely, many-scale) nature of the equations 
of fluctuating hydrodynamics poses many mathematical and computational challenges that are yet 
to be addressed. Notably, it is necessary to develop temporal integrators that can accurately and 
robustly handle the large separation of time scales between different physical processes, such as 
mass and momentum diffusion. The computational techniques we develop here form the foundation 
for incorporating additional physics, such as heat transfer and internal energy fluctuations, phase 
separation and interfacial dynamics, and chemical reactions. 
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We begin Section|Tl]by formulating the fluctuating low Mach number equations for an isothermal 
binary fluid mixture. We present both a traditional pressure (constrained) formulation and a gauge 
(unconstrained) formulation. We analyze the spatio-temporal spectrum of the thermal fluctuations 
in the linearized equations and demonstrate that the low Mach equations eliminate the fast (sonic) 
pressure fluctuations but maintain the correct spectrum of the slow (diffusive) fluctuations. In 



Section III we describe a spatial discretization of the equations that strictly maintains the equation 

we 



of state constraint and also obeys a fluctuation-dissipation balance principle [41] . In Section IV 



develop projected Runge-Kutta schemes for solving the spatially-discretized equations, including a 
midpoint and a trapezoidal second-order predictor-corrector scheme, and a third-order three-stage 
scheme. In Section ?? we study the steady-state spectrum of giant concentration fluctuations 
in the presence of an applied concentration gradient in a mixture of two dissimilar fluids, and 
test the applicability of common approximations that neglect spatial inhomogeneities. In Section 



VI we study the dynamical evolution of giant interface fluctuations during diffusive mixing of 



two dissimilar fluids, using both hard-disk molecular dynamics and low Mach number fluctuating 
hydrodynamics. We find excellent agreement between the two, providing a strong support for the 
usefulness of the fluctuating low Mach number equations as a coarse-grained model of complex 



fluid mixtures. In Section VII we offer some concluding remarks and point out several outstanding 
challenges for the future. 



II. LOW MACH NUMBER EQUATIONS 

The compressible equations of fluctuating hydrodynamics were proposed some time ago jl5j 
and have since been studied and applied successfully to a variety of situations |16j . The presence 
of rapid pressure fluctuations due to the propagation of sound waves leads to stiffness that makes 
it computationally expensive to solve the fully compressible equations numerically, especially for 
typical liquids. It is therefore important to develop fluctuating hydrodynamics equations that 
capture the essential physics in cases where acoustics can be neglected. Developing coarse-grained 
models that only resolve the relevant spatio-temporal scales is a well-studied but still ad hoc 
procedure that requires substantial a priori physical insight [T7] . More precise mathematical mode- 
elimination procedures [421 B3] are technically involved and often purely formal, especially in the 
context of stochastic partial differential equations (SPDEs). 

Here we follow a heuristic approach to constructing fluctuating low Mach number equations, 
starting from the well-known deterministic low Mach equations (which can be obtained via asymp- 
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totic analysis \23\ I24j ) and then adding fluctuations in a manner consistent with fluctuation- 
dissipation balance. We validate the resulting equations by a linearized analysis and comparison to 



the full compressible equations, and in Section VI we validate the low Mach number model through 
a direct comparison to molecular dynamics simulations of diffusive mixing in two dimensions. The 
excellent agreement that we observe demonstrates the utility of the low Mach number equations as 
a coarse-grained model for out-of-equilibrium flows of complex fluid mixtures, at a small fraction 
of the computational cost of particle methods or full compressible simulations. 



A. Compressible Equations 

The starting point of our investigations is the system of isothermal compressible equations of 
fluctuating hydrodynamics for the density p, velocity v, and mass concentration c for a mixture of 
two fluids [16| [T71 \7T\ . In terms of mass and momentum densities the equations can be written as 
local conservation laws, 

dtP + V- (pv) =0 
dt [pv) + V • [pvv'^) = - VP + V • [77 [Vv + V^v) + S] + 

d,{p,) + V- {p,v) =V • [px (Vc + KpVP) + *] , (1) 
where pi — pc is the density of the first component, p2 = (1 - c)p is the density of the second 
component, P{p,c;T) is the equation of state for the pressure at the reference temperature T = 
To = const., g is the gravitational acceleration, and we have neglected bulk viscosity as a weak 
effect. Temperature fluctuations are neglected in this study but can be accounted for with a very 
similar approach. The shear viscosity ry, mass diffusion coefficient X) ^-nd baro-diffusion coefficient 
Kp, can, in general, depend on the state. The baro-diffusion coefficient Kp above [denoted with 
kp/P in Ref. [21], see Eq. (A. 17) in that paper] is not a transport coefficient but rather determined 
from thermodynamics |44j . 

^ {^^Ji|^P)^ ^ _ _ Adp/dc)p ^ {dP/dcl^ 

where p, is the chemical potential of the mixture at the reference temperature and pc = {dp/dc)p, 
and Cy — (dP/dp)^ is the isothermal speed of sound. The capital Greek letters denote stochastic 
momentum and mass fluxes that are (in the absence of bulk viscosity effects) formally modeled as 

E = y/vkpT {W + W^) and * = ^2xpp-^kBTW, (2) 
where kp is Boltzmann's constant, and W and W are standard white- noise random Gaussian tensor 
and vector fields with uncorrelated components. 
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B. Low Mach Equations 

At mesoscopic scales, in typical liquids, sound waves are much faster than momentum diffusion 
and can usually be eliminated from the fluid dynamics description. Formally, this corresponds to 
taking the zero Mach number singular limit oo of the system ([T]). The limiting dynamics 

can be obtained by performing an asymptotic expansion in the Mach number |23j . Heuristically, 
to leading order the low Mach number equations can be obtained by making the anzatz that 
the thermodynamic and mechanical components of pressure can be separated, with the dominant 
thermodynamic component being fixed at a reference state, taken here to be a spatially uniform 
state. This leads to the thermodynamic or equation of state (EOS) constraint 

P(p, c) = Pq = const. (3) 
This constraint means that any change in concentration (equivalently, pi) must be accompanied by 
a corresponding change in density, as would be observed in a system at thermodynamic equilibrium 
held at the fixed reference pressure and temperature. For large-scale systems, in the presence of 
gravity a compressible gas would stratify and it would be more appropriate to take a spatially 
varying Pq that is in mechanical equilibrium with gravity, VPq = pg [45j. Here for simplicity 
we consider a mixture of incompressible fluids of small vertical height. This implies that gravity 
causes negligible changes in density, and the variations in density are dominated by variations in 
composition. Note that we do not account for temperature variations in our isothermal model. 

The equation for pi can be written in primitive (non-conservation) form as the concentration 
equation 

Dc 

p—= pDtC = p{dtc + v - Vc) = V ■ F, (4) 

where the non-advective (diffusive and stochastic) fluxes are denoted with 

F = pxVc + 

By differentiating the EOS constraint along a Lagrangian trajectory we obtain 

^^f3p^^[3V-F = d,p + v-Vp = -p\/-v, (5) 
where the solutal expansion coefficient 

1 fdp 



is determined by the specific form of the EOS and in principle depends on concentration. Equation 
([5]) shows that the EOS constraint can be re-written as a constraint on the divergence of velocity, 

pV - v = -I3V ■ F. 

Note that the usual incompressibility constraint is obtained when the density is not affected by 
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changes in concentration, [3 = 0. When /3 ^ changes in composition (concentration) due to 
diffusion cause local expansion and contraction of the fluid and thus a nonzero V ■ v. 

This analysis leads to the isothermal low Mach number equations for a binary mixture of fluids 
in conservation form, 

dt (pv) + V • {pvv'^) = - Vtt + V • [77 {Vv + V^v) + S] + (6) 
d, (pi) + V • {p,v) =V • F (7) 
d, ip2) + V • (p^v) = V F (8) 
s.t. V - v = - {p-^/3) V F. (9) 
The baro-diffusion term does not appear in the low-Mach equation since it is of thermodynamic 
origin and VPq = 0. The gradient of the non-thermodynamic component of the pressure tt (La- 
grange multiplier) appears in the momentum equation as a driving force that ensures the EOS 
constraint ^ is obeyed. By adding the two density equations (TjS) we get the usual continuity 
equation for the total density, 

dtp + V- (pv) = 0. (10) 
Our conservative numerical scheme is based on Eqs. ( |6|7|9|10| . 

1. Model Equation of State 

In general, the EOS constraint ([S]) is a non-linear constraint. In this work we consider a specific 
linear EOS, 

^ + ^ = '^+^^^^1, (11) 

Pi P2 Pi P2 

where pi and p2 are the densities of the pure component fluids (c = 1 and c = 0, respectively), giving 

P = ^-^= . (12) 

P2 Pi CP2 + (1 - C)pi 

The above form of the density dependence on concentration arises if one assumes that the two 
fluids do not change volume upon mixing. This is a reasonable assumption for liquids that are not 



too dissimilar at the molecular level. Surprisingly the EOS (11) is also valid for a mixture of ideal 
gases, since 

P = Pi+P2 = Po= nkgT = {n, + n^) keT = ( ^ + ) fcsT, 

\mi TO2 / 



where m is molecular mass and n = p/m is the number density. This is exactly of the form (11) 
with pi = miPo/ (kBT) = mrii and p2 = nm2. 
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Even if the specific EOS (11) is not a very good approximation over the entire range of con- 



centration < c < 1, (11) may be a very good approximation over the range of concentrations 
of interest if pi and p2 are adjusted accordingly. In this case pi and p2 are not the densities of 
the pure component fluids but rather fitting parameters that approximate the true EOS in the 
range of concentrations of interest. For small variations in concentration around some reference 
concentration c and density p one can approximate /3 ~ p~^ (dp/dc)- by a constant and determine 



appropriate values of pi and p2 from (12) and the EOS (11) evaluated at the reference state. There 
is therefore no loss of generality by assuming the specific form of the EOS. On the other hand, this 
choice will aid significantly in the construction of simple conservative spatial discretizations that 
strictly maintain the EOS without requiring complicated nonlinear iterative corrections. 



2. Boundary Conditions 

Several different types of boundary conditions can be imposed for the low Mach number equa- 
tions, just as for the more familiar incompressible equations. The simplest case is when periodic 
boundary conditions are used for all of the variables. We briefly describe the different types of 
conditions that can be imposed at a physical boundary with normal direction n. 

For the concentration (equivalently, pi), either Neumann (zero mass flux) or Dirichlet (fixed 
concentration) boundary conditions can be imposed. Physically, a Neumann condition corresponds 
to a physical boundary that is impermeable to mass, while Dirichlet conditions correspond to a 
permeable membrane that connects the system to a large reservoir held at a specified concentration. 
In the case of Neumann conditions for concentration, both the normal component of the diffusive 
flux F„ — and the advective flux piv„ = vanish at the boundary, implying that the normal 
component of velocity must vanish, Vn = 0. For Dirichlet conditions on the concentration, however, 
there will, in general, be a nonzero normal diffusive flux F„ through the boundary. This diffusive 
flux for concentration will induce a corresponding mass flux, as required to maintain the equation 
of state near the boundary. From the condition ([9|, we infer the proper boundary condition for 
the normal component of velocity to be 

v^ = -{p-'P)F,,. (13) 
Note that no additional boundary conditions can be specifled for p since its boundary conditions 
follow from those on c via the EOS constraint. 

For the tangential component of velocity v^, we either impose a no-slip condition = 0, or 
a free slip boundary condition in which the tangential component of the normal viscous stress 
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vanishes, 



dr dn 



In the case of zero normal mass flux, w„ = 0, the free slip condition simplifies to a Neumann 
condition for the tangential velocity, dv^/dn — 0. 

C. Gauge Formalism 

For the purposes of analysis and constructing numerical schemes, it is useful to rewrite the 
constrained low Mach number equations as an initial value problem. In the incompressible case, 
V • f = 0, we can write the constrained Navier-Stokes equations as an unconstrained system by 
eliminating the pressure using a projection operator formalism. The constraint V-i; = is a constant 
linear constraint and independent of the state and of time. However, in the low Mach number 
equations the velocity-divergence constraint V • v = ~l3Dtc depends on concentration, and also 
on time when there are additional (stochastic or deterministic) forcing terms in the concentration 
equation. This makes a projection operator approach cumbersome and a gauge formulation |46| 
an attractive alternative. We note that the gauge formulation presented here is simply a formalism 
that helps us understand the fluctuating low Mach equations and, in particular, the associated 
temporal integrators. 



The density p is not an independent variable due to the EOS constraint (11). We can therefore 
eliminate p from consideration and simply take p = p{c) to be a known function of the concentration. 
With this in mind, the low Mach number system ( 6|7|9 ) has the form 



9tm + Vtt = f {c,v,t) 
dtPi = h{c,v,t) 

V-v = S{c,t), (14) 
where m = pv is the momentum density, and /, h and S are given functions. At present, we will 
assume that these functions are smooth functions of time, which is only justified in the presence of 
stochastic forcing terms in a linearized setting. 
Let us introduce a new variable 

rh — pv — m + V-i/), 



where i/; is a gauge variable, which is not uniquely determined; here we choose the gauge for which 
dt^j = IT with initial condition ijj{t = 0) = 0. The appropriate boundary conditions for V' are linked to 
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the boundary conditions on v; we set ip to be periodic if v is periodic, and employ a homogeneous 



Neumann (natural) boundary condition dip/dn — a Dirichlet condition (13) is specified for the 
normal component of the velocity «„. Note that in the spatially-discrete staggered formulation 
that we employ, the homogeneous Neumann condition follows automatically from the boundary 
conditions on velocity used to define the appropriate divergence and gradient operators in the 
interior of the domain. 

If we know m and ipi we can obtain v — (m — V-^), and write the momentum equation as 
an initial value problem 

dtrh = f{c,v,t) 
dtPi = h{c,v,t) 
V-v = Sic,t). (15) 



By virtue of the relation dtip = tt this system is equivalent to (14). The utility of the gauge 
formulation is that in fact, we do not need to know tp in order to determine v. Therefore, the time 
evolution equation for ip does not actually need to be solved, and in particular, the pressure does 
not need to be computed. To see this, we perform a generalized Hodge decomposition and split 
the velocity into two components, 

V = u + VC, 

where V • m = is a divergence- free (solenoidal or vortical) component, and therefore 

V ■v = V\^S{c,t). 

This is a Poisson problem for ( that is well-posed for appropriate boundary conditions on v. 
Specifically, periodic boundary conditions on v imply periodic boundary conditions for u and C- 



At physical boundaries where a Dirichlet condition (13) is specified for the normal component of 
the velocity, we set u„ — and use Neumann conditions for the Poisson solve, dQ/dn = v„. 
Assume that we know rh ^ pv and want to find v or equivalently u from 

V = u + VC + p^'^Vip. 

Using the condition V • m = we obtain a density-weighted Poisson equation for tp, 

V • (p-'V^p) = - V • (i) - VC) = -V • * + S{c, t). 
Let denote the solution operator to the density-dependent Poisson problem, formally, 

= [v-(p-^v)]-\ 

and also define a density-dependent projection operator Vp defined through its action on a vector 
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field w, 

TpW = w- p-^v [i;' (V • w)] . 

This is a well-known momentum-conserving generalization of the constant-density projection op- 
erator Vw = w — V [V^^ (V • w)\ . We can now write 

u^r,{i- VC) = r,i + p-'v [l;'s{c, t)] - VC. 

This gives 

v = u + VC = ns{i), 

where we introduced an affine transformation IZsic, t) that depends on c and t through p and S{c, t), 
and is defined via its action on a vector field w, 

IZs {w)^w- p-'V [L;^ {V-w-S)]. (16) 
Note that application of 7?.s requires only one Poisson solve and does not actually require computing 

C- 

The calculations above show that by using the gauge formulation we can formally write the low 
Mach number equations in the form of an unconstrained initial value problem 

d,rh = f{c,ns{v),t) (17) 

d,p, - h{c,ns{i),t). (18) 

By adopting the gauge formulation, we can use a method of lines approach for spatially-discretizing 
the system ( 17|18 ), and then apply standard Runge-Kutta temporal integrators to the resulting 



system of ordinary (stochastic) differential equations. 

It is important to point out that the actual independent physical variables in the low Mach 
formulation ( iTp^ ) are the vortical (solenoidal) component of velocity u and the concentration 



c. The velocity v = u + V [V~^S'(c, i)] and density p — p{c) are determined from u and c and the 
constraints; hence they can formally be eliminated from the system. This will become evident 
in the linearized analysis we carry out next, which shows that fluctuations in the vortical velocity 
modes are decoupled from the longitudinal fluctuations. We believe that the gauge formulation can 
be used to analyze the nonlinear equations as well, however, we do not carry out such an analysis 
in this paper. 
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D. Linearized Analysis 

As discussed in more depth in Ref. [22], there are fundamental mathematical difficulties with 
the interpretation of the nonlinear equations of fluctuating hydrodynamics due to the roughness of 
the fluctuating fields. It should be remembered, however, that these equations are coarse-grained 
models with the coarse-graining length scale set by the size of the hydrodynamic cells used in 
discretizing the equations |47] . The spatial discretization removes the small length scales from 
the stochastic forcing and regularizes the equations. It is important to point out, however, that 
imposing such a small-scale regularization (smoothing) of the stochastic forcing also requires a 
suitable renormalization of the transport coefficients [IJ |l8j , as we discuss in more detail in Section 

ED 

As long as there are sufficiently many molecules per hydrodynamic cell the fluctuations in the 
spatially-discrete hydrodynamic variables will be small and the behavior of the nonlinear equations 
will closely follow that of the linearized equations of fluctuating hydrodynamics [22], which can be 
given a precise meaning [59]. It is therefore crucial to understand the linearized equations from 
a theoretical perspective, and to analyze the behavior of the numerical schemes in the linearized 
setting [IT] . 

1. Equilibrium Fluctuations: Compressible Equations 

Some of the most important quantities predicted by the fluctuating hydrodynamics equations 
are the equilibrium structure factors (static covariances) of the fluctuating fields. These can be 
obtained by linearizing the compressible equations ([T]) around a uniform reference state, p = po + Sp, 
c — Co + Sc, V — Sv, P — + SP where 

SP^cl [{6p)-Pp{Sc)], 

and then applying a spatial Fourier transform \W\ UT] . Owing to fluctuation-dissipation balance 
the static structure factors are independent of the wavevector k at thermodynamic equilibrium, 

5„,„(fc)= ((H(H*) 

Note that density fluctuations do not vanish even in the incompressible limit Ct ^ oo unless (3 = 0. 
While fluctuations in pi and p2 are uncorrelated, the fluctuations in concentration and density are 



I2 ' P 



Po ^ksTo I 



(19) 
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correlated even at equilibrium, 

Sc,p = {{^p) ('^c) ) = ^^^^7^ ^ Po/3<S'c,c- 
We will see below that the low Mach equations correctly reproduce the static covariances of density 
and concentration in the limit ct — > oo. 

The dynamics of the equilibrium fluctuations can also be studied by applying a Fourier-Laplace 
transform in time in order to obtain the dynamic structure factors (equilibrium correlation func- 
tions) as a function of wavenumber k and wavefrequency w |16| I41| . It is well-known that the 
dynamic spectrum of density fluctuations Sp^p {k,uj) exhibits three peaks for a given k, one central 
Rayleigh peak at small frequencies (slow concentration fluctuations), and two symmetric Brillouin 
peaks centered around lo « crk. As the fluid becomes less compressible (i.e., the speed of sound 
increases), there is an increasing separation of time-scales between the side and central spectral 
peaks. As we will see below, the low Mach equations reproduce the central peaks in the dynamic 
structure factors only, eliminating the side peaks and the associated stiff dynamics. 



2. Low Mach Equations 

We now examine the spatio-temporal correlations of the steady-state fluctuations in the low 
Mach number equations ( |6|7[9|10l ). In order to model the nonequilibrium setting in which giant 
concentration fluctuations are observed, we include a constant background concentration gradient 
in the equations. Note that a density gradient will accompany a concentration gradient, and this 
can introduce some additional terms in F depending on how px depends on concentration. For 
simplicity, we assume px is a constant so that the diffusive term V • F in ([T]) is simply pxV^c. We 
also assume the viscosity r/ is spatially constant, to get the simplified coupled velocity-concentration 
equations, 

DtV = - p'^Vtt + uV'^v + p~^ (V • S) + gr 
Dtc =xV^c + p"' (V • *) 
V ■v=- f3DtC. (20) 



where v = ij/p and p = p{c) is given by ( 11 ) 



We linearize the equations (20) around a steady state, c — c+Sc, v — v + Sv = Sv, and tt = tt + Stt, 
where the reference state is in mechanical equilibrium, p^-^Vn = g. We denote the background 
concentration gradient with h = Vc, and drop the bars assuming the mean concentration (density) 
does not change much across the domain, allowing for a quasi-periodic or weak-gradient approxi- 
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mation \36\ I50j. In the linear approximation, the EOS constraint relates density and concentration 
fluctuations, Sp = pl3{Sc). The term v ■ is second order in the fluctuations and drops out, but 
the advective term v ■ Vc leads to a term {5v) ■ h in the concentration equation. The forcing term 
due to gravity becomes p^^ (Sp) g = /3 {6c) g. After a spatial Fourier transform, the linearized form 



of (20) becomes a collection of stochastic differential equations, one system of linear additive- noise 

equations per wavenumber, 

dt(Sv^ = -ip-'k(^^'^ -ve (Sv^ +ip-'k-'E + Pg(jc)j (21) 
dt (Sc) = h ■ (&v) - xfc' (Sc) + ip-^ (k ■ (22) 



k- (5v] ^ -f3 



ixk (^c) + p-' (fc • . (23) 
Replacing the right hand side of (23) with zero leads to the incompressible approximation used in 
Ref . [36J , corresponding to the Boussinesq approximation of taking the limit (3-^0 while keeping 
the product f5g constant. 

3. Equilibrium Fluctuations: Low Mach Equations 

Let us first compare the dynamics of the equilibrium fluctuations (h — O) in the low Mach 
equations with those in the complete compressible equations. For simplicity of notation we will 
continue to use the hat symbol to denote the space-time Fourier transform. 

In the wavenumber-frequency (fe,w) Fourier domain, the concentration fluctuations in the ab- 



sence of a gradient are obtained from (22) 



iuj + xk V 

which is the same as the compressible equations. The density fluctuations follow the concentration 
fluctuations, 5p = pPSc, and the dynamic structure factor for density shows the same central 
Rayleigh peak as obtained from the isothermal compressible equations |16j . 



where we used Eq. (^2| for the covariance of ^ . This shows that the low Mach number equations 



correctly reproduce the slow fluctuations (small w) in density and concentration, while eliminating 
the side Brillouin peaks associated with the fast isentropic pressure fluctuations. 

The fluctuations in velocity, however, are different between the compressible and low Mach 
number equations. Let us first examine the transverse (solenoidal) component of velocity Sv^ — V5v, 
where V is the constant-density orthogonal projection onto the space of divergence-free velocity 
fields, V — I — k~'^{kk*) in Fourier space. Applying the projection operator to the velocity equation 
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(21) shows that the fluctuations of the solenoidal modes are the same as in the incompressible 



approximation , 



dt (Sv, ] ^~i,k^( 6v, ) + ip-^k ■ T'S + l3Vg ( 5c 



The fluctuations of the compressive velocity component dvi ~ k - (Sv^ , on the other hand, are driven 
by the stochastic mass flux '4', as seen from eq. (23) at thermodynamic equilibrium, 

' iu! + ^ " " 

The dynamic structure factor (space-time Fourier spectrum) of the longitudinal component 

does not decay to zero as a; -h- oo. This indicates that the fluctuations of velocity are not only 
white in space but also white in time. In the incompressible approximation /3 -> the longitudinal 
velocity fluctuations vanish and the static spectrum of the velocity fluctuations is equal to the 
projection operator, „ = V [22j. In the compressible equations, the dynamic structure factor for 
the longitudinal component of velocity decays to zero as w — > oo because it has two sound (Brillouin) 
peaks centered around w w cxk, in addition to the central diffusive (Rayleigh) peak. The low Mach 
number equations reproduce the central peak (slow fluctuations) correctly, replacing the side peaks 
with a flat spectrum for large w. The origin of this unphysical behavior is the unjustifled assumption 
of infinite separation of time scales between the propagation of sound and the diffusion of mass, 
momentum and energy. In reality, the same molecular motion underlies all of these processes and 
the incompressible or the low Mach number equations cannot be expected to reproduce the correct 
physical behavior at very short time scales (w > crk). 

The fact that the velocity fluctuations are white in space and in time poses a further challenge 
in interpreting the nonlinear low Mach number equations, and in particular, numerical schemes 
may not converge to a sensible limit as the time step goes to zero. In practice, just as the spatial 
discretization of the equations imposes a spatial smoothing or regularization of the fluctuations, the 
temporal discretization of the equations imposes a temporal smoothing and filters the problematic 
large frequencies. In the types of problems we study in this work the equilibrium concentration 
fiuctuations can be neglected, '4' « 0, because the concentration fluctuations are dominated by 
nonequilibrium effects. This eliminates the problematic white component in the velocity fiuctu- 
ations. In the more general setting we believe the gauge formulation can be used to write well- 
behaved stochastic differential equations for the vortical component of the velocity u = v — VC and 
concentration c, and also analyze the accuracy with which temporal integrators approximate those 
equations. 
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4- Nonequilibrium Fluctuations: Low Mach Equations 



If we neglect the term involving "if in (23) and eliminate the Lagrange multiplier (non- 



thermodynamic pressure) n using (23), we obtain the linearized velocity equation in Fourier space 



ip-^k-VS + 13 [Scj Vg 
h ■ (Ju) k + i/3x {v - X) k'^ (Jc) k. (24) 
It is straightforward to obtain the steady-state covariances (static structure factors) in the pres- 
ence of a concentration gradient from the linearized system of velocity-concentration equations 
( 22|24 ) [41J. The procedure amounts to solving a linear system for three covariances (velocity- 
velocity, concentration-concentration, and velocity-concentration). These types of calculations are 
particularly well-suited for modern computer algebra systems like Maple and can be carried out 
for arbitrary wavenumber and background concentration gradient. We omit the full solution for 
brevity. 

Experiments measure the steady-state spectrum of concentration fluctuations averaged along 
the gradient [21 HO], and we will therefore focus on wavenumbers perpendicular to the gradient, 
k ■ h — 0. A straightforward calculation shows that the concentration fluctuations are enhanced as 
the square of the applied gradient. 



el component relative to gravity, respectively. 



(25) 



where _L and || denote the perpendicular and parall 
The term in the denominator involving comes from the low Mach number constraint (|9| and 
is usually negligible since the concentration gradient is parallel to gravity or xl^ ^ 1- Without 



this term the result (25) is the same result as obtained in [36], and shows that fluctuations at 



wavenumbers below fcj^ — (i^x) are suppressed by gravity, as we study numerically in Section 



?? 



III. SPATIAL DISCRETIZATION 



Our spatio-temporal discretization follows a "method of lines" approach in which we first dis- 
cretize the equations ( |6|7|9|10l ) in space and then integrate the resulting semi-continuum equations 
in time. The spatial discretization we employ follows closely the spatial discretization of the 
constant-coefficient incompressible equations described in Ref. |22j . Therefore, we focus here on 
the differences, specifically, the use of conserved variables, the handling of the variable-density pro- 
jection and variable-coefficient diffusion, and the imposition of the low Mach number constraint. 
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Figure 1: Staggered (MAC) finite- volume discretization on a uniform Cartesian two-dimensional grid. {Left) 
Control volume and flux discretization for cell-centered scalar fields, such as densities p and pi. (Middle) 
Control volume for the a;-component of face-centered vector fields, such as (Right) Control volume for 
the j/-component of face-centered vector fields, such as my. 

Note that the handling of the stochastic momentum and mass fluxes is identical to that described 



For simplicity of notation, we focus on two dimensional problems, with straightforward gener- 
ahzation to three spatial dimensions. Our spatial discretization follows the commonly- used MAC 
approach in which the scalar conserved quantities p and pi are defined on a regular Cartesian 
grid. The vector conserved variables m — pv are defined on a staggered grid, such that the fc"" 
component of momentum is defined on the faces of the scalar variable Cartesian grid in the fc**" di- 
rection, see Fig. [l] For simplicity of notation, we often denote the different components of velocity 
as V = {u,v) in two dimensions and v — {u,v,w) in three dimensions. The terms "cell-centered", 
"edge-centered", and "face-centered" refer to spatial locations relative to the underlying scalar grid. 
Our discretization is based on calculating fluxes on the faces of a finite-volume grid and is thus 
locally conservative. It is important to note, however, that for the MAC grid different control 
volumes are used for the scalars and the components of the momentum, see Fig. [T] 

From the cell-centered p and pi we can define other cell-centered scalar quantities, notably, the 
concentration c^^ = iPi),j/Pi,j and the transport quantities Xi,] and riij, which typically depend 
on the local density pij and concentration dj (and temperature for non-isothermal models), and 
can, in general, also depend on the spatial position of the cell {x,y) = (iAx, jAy). In order to 
define velocities we need to interpret the continuum relationship m = pv on the staggered grid. 
This is done by defining face-centered scalar quantities obtained as an arithmetic average of the 
corresponding cell-centered quantities in the two neighboring cells. Specifically, we define 



inRef. 



Pi,j + Pi+l,j 

2 



u. 




(26) 
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except at physical boundaries, where the value is obtained from the imposed boundary conditions 



(see Section III B 3 ) . Arithmetic averaging is only one possible interpolation from cells to faces 
|52| . In general, other forms of averaging such as a harmonic or geometric average or higher- 
order, wider stencils [4H [53] can be used. Most components of the spatial discretization can easily 
be generalized to other choices of interpolation. As we explain later, the use of linear averaging 
simplifies the construction of conservative advection. 



A. Diffusion 



In this section we describe the spatial discretization of the diffusive mass flux term V • px^c in 
Q . The discretization is based on conservative centered differencing |4H IM] , 



(V • pxVc).,, = Ax-i 
where, for example. 



dc 



PX 



dc\ 



dx ) ,1, . dx ) . I, . 



PX 



dc 
dy 



J+V2 



PX 



dc 
dy 



PX 



dc 
dx 



(A+V2,j) (X» + V2j) 



(27) 



(28) 



and Xi+Va.j interpolated face-centered diffusion coefficient, for example, as done for p in Eq. 



(26) 



Xi.j Xi + l,J 

Xi+'^h.i ^ 2 ' 

except at physical boundaries, where the value is obtained from the imposed boundary conditions. 

Regardless of the specific form of the interpolation operator, the same face-centered diffusion 
coefficient Xi+Va.j ™ust be used when calculating the magnitude of the stochastic mass flux on face 



This ensures discrete fluctuation-dissipation balance in the linearized setting. Specifically, at ther- 
modynamic equilibrium the static covariance of the concentration is determined from the equi- 
librium value of (pPc^) (thermodynamics) independently of the particular values of the transport 



coefficients (dynamics), as seen in (19) and dictated by statistical mechanics principles. 



B. Viscous Terms 



In Ref. a Laplacian form of the viscous term -qV^v is assumed, which is not applicable when 
viscosity is spatially varying and V ■ v — S ^ 0. In two dimensions, the divergence of the viscous 
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stress tensor in the momentum equation neglecting bulk viscosity effects, is 



2-^ 

dy 



(29) 



The discretization of the viscous terms requires 77 at cell-centers and edges (note that in two 
dimensions the edges are the same as the nodes (i + 1/2, j + 1/2) of the grid). The value of 77 at a 
node is interpolated as the arithmetic average of the four neighboring cell-centers, 

except at physical boundaries, where the values are obtained from the prescribed boundary condi- 
tions. The different viscous friction terms are discretized by straightforward centered differences. 



Explicitly, for the x-component of momentum 



d_ 
dx 



du 
^ dx 



Ax-' 



du 
^ dx 



i + l.j 



du 



with 



du 
^dx 



Ax 



Similarly, for the term involving a second derivative in y, 



d ( du\ 



dy \^ dy J 



^Ay-' 



i + V2.j 



du 
dy 



i + V2,J + V2 



du 
dy 



i + V2,i-V2. 



with 



du\ 



+ V2,i + V2 



'?i+V2,J + V2 



Ay 



A similar construction is used for the mixed-derivative term. 



d 



dv\ 



dy \ dx J 



^Ay-' 



i + V2.j 



dv\ 
^d-x) 



i + V2,i + V2 



dv\ 
dx I 



i + V2,J-V2. 



with 



dv 



^i+V2J + V2 



i + l.i+72 



■ V. 



Aa 



i + V2,i + V2 

The stochastic stress tensor discretization is described in more detail in Ref. [22J and applies 
in the present context as well. For the low Mach number equations, just as for the compressible 
equations, the symmetric form of the stochastic stress tensor must be used in order to ensure discrete 
fluctuation-dissipation balance between the viscous dissipation and stochastic forcing. Additionally, 
when Tj is not spatially uniform the same interpolated viscosity ryj,,.!,^ as used in the viscous terms 
must be used when calculating the amplitude in the stochastic forcing y/rjk^ at the edges (nodes) 
of the grid. 
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Advection 



It is challenging to construct spatio-temporal discretizations that conserve the total mass while 
remaining consistent with the equation of state [2S\ [501 ES] j as ensured in the continuum context 



by the constraint (|9|) . We demonstrate here how the special linear form of the constraint (11) can 



be exploited in the discrete context. Following Ref. [22], we spatially discretize the advective terms 
in ([7|) using a centered (skew-adjoint ^56j) discretization, 



+ V2,i ■"' + V2,i 



-V2,J 



i-V2 ^«.j-y2 



(30) 



and similarly for (10). We would like this discrete advection to maintain the equation of state ( |11[ ) 
at the discrete level, that is, maintain the constraint relating (/Oi);^ and (^2)^^ in every cell 

Because the different dimensions are decoupled and the divergence is simply the sum of the 
one-dimensional difference operators, it is sufficient to consider ([7]) in one spatial dimension. The 
method of lines discretization is given by the system of ODEs, one differential equation per cell i, 



(Pl), + 1/,U>+V2 - iPl)^- 



V2"»-V2 



and similarly for {dtP2)-. As a shorthand, denote the quantity that appears in (11) with 

Pi P2 



1. 



If we use the linear interpolation (26) to calculate face-centered densities, then because of the 



linearity of the EOS the face-centered densities obey the EOS if the cell-centered ones do, since 
The rate of change of S in cell i is 



= 1. 



-y2 - P^-y2) 



Axid,S\ = (p-^/3) (F,+v, 

= {p''/3) (F,+ y, - - {u,^y^ - = 0. 

This simple calculation shows that the EOS constraint J = 1 is obeyed discretely in each cell at 
all times if it is initially satisfied and the velocities used to advect mass obey the discrete version 
of the constraint ([O]), 



' (^«+V2,. - "«-V2,.) + ' (^..j+V2 - ^^..-Vi) (31) 
1 1 

.Pl P2, 

in two dimensions. Our algorithm ensures that advective terms are always evaluated using a 



V2,jJ 



Ay-' {F^ 



3-\ 



)] 



discrete velocity field that obeys this constraint. This is accomplished by using a discrete projection 
operator, as we describe in the next section. 

The spatial discretization of the advection terms in the momentum equation ([6]) is constructed 
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using centered differences on the corresponding shifted (staggered) grid, as described in Ref. 
For example, for the x-component of momentum m^, = pu, 

where simple averaging is used to interpolate momenta to the cell centers and edges (nodes) of the 
grid, for example, 

M..^ = (^.).„..„ = ( ^"^-^'--^-;^""^'--^-- ) 1"'-^"- ) . (33) 

Because of the linearity of the interpolation procedure, the interpolated discrete velocity used to 
advect obeys the constraint (31) on the shifted grid, with a right-hand side Si^i^^ j interpolated 



using the same arithmetic average used to interpolate the velocities. In particular, in the incom- 
pressible case all variables, including momentum, are advected using a discretely divergence-free 
velocity, ensuring discrete fluctuation-dissipation balance [22l [5l] . 

It is well-known that the centered discretization of advection we employ here is not robust 
for advection-dominated flows, and higher-order limiters and upwinding schemes are generally 
preferred in the deterministic setting |57j . However, these more robust advection schemes add 
artificial dissipation, they lead to a violation of fluctuation-dissipation balance. In Appendix |^ 
we describe an alternative filtering procedure that can be used to handle strong advection while 
continuing to use centered differencing. 



2. Discrete Projection 



We now briefly discuss the spatial discretization of the affine operator TZs defined by (16), as 
used in our explicit temporal integrators. The discrete projection takes a cell-centered discrete 
velocity field v = {u, v) and a velocity divergence S and projects v = TZg {v) onto the constraint 



(31 ) in a conservative manner. Specifically, the projection consists of finding a cell-centered discrete 



scalar field such that 

pv = pv — V(j), and V ■ v = S, 
where the gradient is discretized using centered differences, e.g.. 
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The pressure correction </> is the solution to the variable-coefficient discrete Poisson equation, 



1 
Ax 
1 

'Ay 



1 



1 



Ax 



1 



Ay 



1 

A,,-V2 



~ Vi.j-l 

Ax 

i>i,j ~ (f'i.j-1 



Ax 



Ay 

- V 

a7 



(35) 



which can be solved using a standard multigrid approach \52\ . 



3. Boundary Conditions 

The handling of different types of boundary conditions is relatively straightforward when a 
staggered grid is used and the physical boundaries are aligned with the cell boundaries for the 
scalar grid. Interpolation is not used to obtain values for faces, nodes or edges of the grid that lie 
on a physical boundary, since this would require "ghost" values at cell centers lying outside of the 
physical domain. Instead, whenever a value of a physical variable is required at a face, node, or 
edge lying on a physical boundary, the boundary condition is used to obtain that value. Similarly, 
centered differences for the diffusive and viscous fluxes that require values outside of the physical 
domain are replaced by one-sided differences that only use values from the interior cell bordering 
the boundary and boundary values. 

For example, if the concentration is specified at the face (i + j), the diffusive flux discretization 



( 28 ) is replaced with 



-iV ■ V Ax/2 

where c.+y^j is the specified boundary value, the density Pi+y^.j is obtained from c,^!/^ ^ using the EOS 
constraint, and the diffusion coefficient Xi+'^k.j is calculated at the specified values of concentration 
and density. Similar straightforward one-sided differencing is used for the viscous fluxes. As 
discussed in Ref. [22j , the use of second-order one-sided differencing is not required to achieve global 
second-order accuracy, and would make the handling of the stochastic fluxes more complicated 
because it leads to a non-symmetric discrete Laplacian. Note that for the nonlinear low Mach 
number equations our approach is subtly different from linearly extrapolating the value in the 
ghost cell Ci+ij = 2cj_,.y2 J — c^. Namely, the extrapolated value might be unphysical, and it might 
not be possible to evaluate the EOS or transport coefficients at the extrapolated concentration. For 
Neumann-type or zero-flux boundary conditions, the corresponding diffusive flux is set to zero for 
any faces of the corresponding control volume that lie on physical boundaries, and values in cells 
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outside of the physical domain are never required. The corresponding handhng of the stochastic 
fluxes is discussed in detail in Ref. |22j . 

The evaluation of advective fluxes for the scalars requires normal components of the velocity at 
the boundary. For faces of the grid that lie on a physical boundary, the normal component of the 



velocity is determined from the value of the diffusive mass flux at that face using (13). Therefore, 
these velocities are not independent variables and are not solved for or modified by the projection 
71s ■ Specifically, the discrete pressure (/) is only defined at the cell centers in the interior of the grid. 



and the discrete Poisson equation (35) is only imposed on the interior faces of the grid. Therefore, 
no explicit boundary conditions for cj) are required when the staggered grid is used, and the natural 
homogeneous Neumann conditions are implied. Advective momentum fluxes are only evaluated on 
the interior faces and thus do not use any values outside of the physical domain. 

IV. TEMPORAL INTEGRATION 



Our temporal integrators are based on the gauge formulation ( 17|18 ) of the low Mach equations. 



The gauge formulation is unconstrained and enables us to use standard temporal integrators for 

initial-value problems. In the majority of this section, we assume that all of the fields and differential 

operators have already been spatially discretized and focus on the temporal integration of the 

resulting initial-value problem. 

Because in the present schemes we handle both diffusive and advective fluxes explicitly, the 

time step size At is restricted by well-known CFL conditions. For fluctuating hydrodynamics 

applications the time step is typically limited by momentum diffusion, 

_ vAt 1 
"■^ " ^ 2d' 

where d is the number of spatial dimensions. The design and implementation of numerical methods 
that handle momentum diffusion semi-implicitly, as done in Ref. |22j for incompressible flow, is 
substantially more difficult for the low Mach number equations and will be the subject of future 
work. 

Our temporal discretization will make use of the special form of the EOS and the discretization 



of mass advection described in Section IIIB 1 in order to strictly maintain the EOS relation (11) 
between density and concentration in each cell at all intermediate values. Therefore, no additional 
action is needed to enforce the EOS constraint after an update of pi and p. This is, however, only 
true to within the accuracy of the Poisson solver and also roundoff, and it is possible for a slow 



drifting off the EOS to occur over many time steps. In Section IV C, we describe a correction 



24 



that prevents such drifting and ensures that the EOS is obeyed at all times to essentially roundoff 
tolerance. For simplicity, we will sometimes omit the explicit update for the density p and instead 
focus on updating pi and the momentum density m — pv, with the understanding that p is updated 
whenever pi is. 



A. Euler Scheme 



The foundation for our higher-order explicit temporal integrators is the first-order Euler method 
applied to the gauge formulation ( |17|18| ) . The specific forms of the functions appearing in ( 17|18 ) 
for the low Mach equations ( |6|7|9|10l ) are given by 

h = V ■ i-piV + F) 

/ = V • [-pvv"^ + {Vv + V^t>) pg 
S = - {p-'p) V • F 
F = pxVc + *, 



where all of the terms have already been spatially discretized as described in Section III 



1. Gauge-Free Euler Update 

We use a superscript to denote the time step and the point in time where a given term is 
evaluated, e.g., /i" = h{c",v",t"), and we denote the time step size with At = — t". Assume 
that at the beginning of timestep n we know rh" and define 



to be the velocity at time nAt obtained by enforcing the constraint (14). Here T^.^ denotes the affine 



transformation ( 16 ) with all terms evaluated at the beginning of the time step, so that V • = S". 
An Euler step for the low Mach equations then consists of the update 

Tn"+i = m" + At/", (36) 
together with an update of the density p"^^ consistent with p"'^^. 

At the beginning of the next time step, will be calculated from m"'^^ by applying TZg'^^, 
and it is only ?;"+-^ that will actually be used during time step n + 1. We therefore do not need to 
explicitly store m"'^^ and can instead replace it with Tn"+^ = without changing any of the 

observable results. This is related to the fact that the gauge is de facto arbitrary and, in the present 
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setting, the gauge formulation is simply a formalism to put the equations in an unconstrained form, 
without any physical significance. The difference between rh and m is a (discrete) gradient of a 
scalar. Since our temporal integrators only use linear combinations of the intermediate values, the 
difference between the final result for rh"'^^ and m" is also a gradient of a scalar and replacing 
m"+^ with m"+^ simply amounts to redefining the (arbitrary) gauge variable. For these reasons, 
we use the following Euler advance as the foundation of our temporal integrators, 

(37) 



(p"+^)"'(m" + At/") 



which is equivalent to (36) with proper initialization. 



2. Stochastic Forcing 



Thermal fluctuations cannot be straightforwardly incorporated in (37) because it is not clear 



how to define 'JZg'^^. In the deterministic setting, S* is a function of concentration and density and 
can be evaluated pointwise at time level n + 1. When the white- in-time stochastic concentration 
flux * is included, however, S cannot be evaluated at a particular point of time. Instead, one must 
think of as representing the average stochastic flux over a given time interval St, which can be 
expressed in terms of the increments \/6iW of the underlying Wiener processes, 

* (^St, = 



StAV 



where W is a collection of normal variates generated using a pseudo-random number generator, 
and AV^ is the volume of the hydrodynamic cells. Similarly, the average stochastic momentum flux 
over a time step is modeled as 

where W are normal random variates. As described in more detail in Ref. [22], stochastic fluxes 
are spatially discretized by generating normal variates on the faces of the grid on which the corre- 
sponding variable is discretized, independently at each time step. 

With this in mind, we first evaluate the velocity divergence associated with the constraint using 
the particular sample of '4', 



S=- {p-^13) V • pxVc + * Vf) 



We then define a discrete affine operator TZp (St, w) in terms of its action on the momentum m 



IZf (st, w'j (m) = pTLs (p^'m) . 



Using this shorthand notation, the momentum update in (37) in the presence of thermal fluctuations 
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can be written as 



7?.^+' ( At W 



(m" + Atr). 



Observe that this is a conservative momentum update since the apphcation of TZp subtracts the 
(discrete) gradient of a scalar from the momentum. In actual implementation, it is preferable to 
apply 'R-p'^^ at the beginning of the time step n + 1 instead of at the end of time step n, once the 
value 5*"+^ is computed from the diffusive and stochastic fluxes for the concentration. 

3. Euler-Maruyama Update 



m' = 



Following the above discussion, we can write an Euler-Maruyama temporal integrator for the 
low Mach number equations in the shorthand notation, 

p'l + At h'' + (Ai, W") 
m"+^ = m" + A^ /" + /" {At, VF") , (38) 
where W" and W are collections of standard normal variates generated using a pseudo-random 
number generator independently at each time step. Here the deterministic increments are written 
using the shorthand notation, 

7 = V • [-pvv"^ + rj [Vv + V^t>)] + pg 
h = V ■ {-piV + px'^c) . 
The stochastic increments are written in terms of 

/ {St, W) = [V • S {St, W)] St = V 



V {ksT) St 



AV 



{W + W'' 



(st, w) 



St = V 



2xPKHkBT)St 



AV 



W 



where W and W are Gaussian variables whose moments match expectation values of stochastic 
integrals of the underlying Wiener processes to certain order in At |54j . 

Before describing higher-order temporal integrators, we summarize our implementation of a 



single Euler step (38). This forms the core procedure which the higher-order Runge-Kutta schemes 



call several times during one time step. 

1. Generate the vectors of standard Gaussian variates W" and W 



2. Calculate diffusive and stochastic fluxes for pi using (28), 

= (pxVc)" + (a*, W 
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3. Solve the Poisson problem (35) with 



5" = - ( ^ - — 1 V -F" 

.Pi P2, 



to obtain the velocity v" from v" = rh^/p" using (34), enforcing V • = S". 



4. Calculate viscous and stochastic momentum fluxes using (29), 



V • [rj (Vv + V^v)Y + V • [S" (At, W")] . 

5. Calculate external forcing terms for the momentum equation, such as the contribution —p"g 
due to gravity. 



6. Calculate advective fluxes for mass and momentum using (30) and (32). 



7. Update mass and momentum densities, including advective, diffusive, stochastic and external 
forcing terms, to obtain /o"^^ and m"+^. Note that this update preserves the EOS 



constraint as explained in Section III B 1 



B. Higher- Order Temporal Integrators 



A good strategy for composing higher-order temporal integrators for the low Mach number 



equations is to use a linear combination of several projected Euler steps of the form (38). In this 



way, the higher-order integrators inherit the properties of the Euler step. In our case, this will be 
very useful in constructing conservative discretizations that strictly maintain the EOS constraint 
and only evaluate fluxes at states that strictly obey the EOS constraint. 

The incorporation of stochastic forcing in the Runge-Kutta temporal integrators that we use 
is described in Refs. |4H I54|: here we only summarize the resulting schemes. We note that the 
stochastic terms should be considered additive noise, even though we evaluate them using an 
instantaneous state like multiplicative noise |22j . 



1. Explicit Midpoint Rule 



A weakly second-order temporal integrator for ( 17|18 ) is provided by the explicit midpoint rule, 
which can be summarized as follows. First we take a projected Euler step to estimate midpoint 
values (denoted here with superscript n + 1/2) , 
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m 



Pi 



1+V2 



= m 



At 
At 

T 

At 



h" + h" 

r + / 



At 



At 



(39) 



and then we complete the time step with another Euler-hke update 



7^^+'''^ (At, w 

Ath*'"+'^^ + h 
-At~f'"'''^" + 



■"+'''^ (At, W 



(At, W") , 



(40) 



where the standard Gaussian variates 



V2 



and the vectors of standard normal variates and W2 are independent, and similarly for W" 
and Note that and are used in 6ot/i the predictor and the corrector stages, while Wj 
and are used in the corrector only. Physically, the random numbers I\f2 (and similarly for 



W ^ ) correspond to the increments of the underlying Wiener processes ABi = At/2 W1 over the 
first half of the time step, and the random numbers W^/^ correspond to the Wiener increments 
for the second half of the timestep |54j . 

Note that both the midpoint and the endpoint values for density and concentration obey the 
EOS. We numerically observe that the midpoint rule does not exhibit a systematic numerical drift 



in the EOS, and can therefore be used without the correction procedure described in Section IV C 



The analysis in Ref. [5jJ indicates that for the incompressible case the midpoint scheme exhibits 
second-order weak accuracy in the nonlinear setting. Furthermore, in the linearized setting it 
reproduces the steady-state covariances of the fluctuating fields to third order in the time step size. 



2. Explicit Trapezoidal Rule 

An alternative second-order scheme is the explicit trapezoidal rule, in which we first take a 
predictor Euler step 
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The corrector step is a linear combination of the predictor and another Euler update, 
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(43) 
(44) 



At/i*'"+^ + (At, W 

m-' ' ^ + At Z*'"^' + (At, T^")" 

and reuses the same random increments as the predictor step. 

Note that both the predicted and the corrected values for density and concentration obey the 
EOS. We numerically observe that the trapezoidal rule does exhibit a slow but systematic numerical 
drift in the EOS, and therefore it is necessary to use the correction procedure described in Section 



IV C at the end of each time step. The analysis in Ref. [54J indicates that for the incompressible 



case the trapezoidal scheme exhibits second-order weak accuracy in the nonlinear and linearized 
settings. 



3. Three- Stage Runge-Kutta (RK3) Rule 



We have also tested and implemented the three-stage Runge Kutta scheme that was used in 
Refs. \12\ I41j . This scheme can be expressed as a linear combination of three Euler steps. The 
first stage is a predictor Euler step, 

(At, W")] (m") 

At/i" + fAt, W^"") (45) 

(46) 
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The second stage is a midpoint predictor 
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and a final corrector stage completes the time step 
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Here the stochastic fluxes between difl'erent stages are related to each other via 

(2V2 + V5) 

W" =W" + ^ -^W" 

5 

(-4^2 + 3 V3) 
5 



W**'''=W'^ + ^ — ^^^^ 
where T^" and are independent and generated independently at each RK3 step, and similarly 
for W. The weights of W2 are chosen to maximize the weak order of accuracy of the scheme while 
still using only two random samples of the stochastic fluxes per time step |54] . 

The RK3 method is third-order accurate deterministically, and stable even in the absence of 
diffusion/viscosity (i.e., for advection-dominated flows). Note that the predicted, the midpoint 
and the endpoint values for density and concentration all obey the EOS. We numerically observe 
that the RK3 scheme does exhibit a systematic numerical drift in the EOS, and therefore it is 



necessary to use the correction procedure described in Section IV C at the end of each time step. 
The analysis in Ref. |54j indicates that for the incompressible case the RK3 scheme exhibits 
second-order weak accuracy in the nonlinear setting. In the linearized setting it reproduces the 
steady-state covariances of the fluctuating fields to third order in the time step size. 



C. EOS drift 



While in principle our temporal integrators should strictly maintain the EOS, roundoff errors 
and the finite tolerance employed in the iterative Poisson solver lead to a small drift in the constraint 
that can, depending on the specific scheme, lead to an exponentially increasing violation of the EOS 
over many time steps. In order to maintain the EOS at all times to within roundoff tolerance, we 
periodically apply a globally-conservative L2 projection of p and pi onto the linear EOS constraint. 

This projection step consists of correcting pi in cell k using 

fc' 

where N is the number of hydrodynamic cells in the system and 

pI + pV pI + pI 

Note that the above update, while nonlocal in nature, conserves the total mass Yl,y (Pi)fc'- ^ similar 
update applies to p2 , or equivalently, p — pi + P2- 
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V. GIANT CONCENTRATION FLUCTUATIONS 

Advection of concentration by thermal velocity fluctuations in the presence of large concentra- 
tion gradients leads to the appearance of giant fluctuations of concentration, as has been studied 
theoretically and experimentally for more than a decade [21 [36l HOl I58j . These giant fluctuations 
were previously simulated in the absence of gravity in three dimensions by some of us in Ref . [22j , 
and good agreement was found with experimental results [2j. In those previous studies the incom- 
pressible equations were used, that is, it was assumed that concentration was a passively-advected 
scalar. However, it is more physically realistic to account for the fact that the properties of the fluid, 
notably the density and the transport coefficients, depend on the concentration. In Ref. [40J a series 
of experiments were performed to study the temporal evolution of giant concentration fluctuations 
during the diffusive mixing of water and glycerol, starting with a glycerol mass fraction of c = 0.39 
in the bottom half of the experimental domain, and c = in the top half. Because it is essentially 
impossible to analytically solve the full system of fluctuating equations in the presence of spatial 
inhomogeneity and nontrivial boundary conditions, the existing theoretical analysis of the diffusive 
mixing process [36] makes a quasi-periodic constant-coefficient incompressible approximation. 

For simplicity, in this section we focus on a time-independent problem and study the spectrum 
of steady-state concentration fluctuations in a mixture under gravity in the presence of a constant 
concentration gradient. This extends the study reported in Ref. |22J to account for the fact that 
the density, viscosity, and diffusion coefficient depend on the concentration. For simplicity, we 
do two-dimensional simulations, since there are no qualitative differences between the spectrum 
of concentration fluctuations in two and three dimensions [22] (note, however, that in real space, 
unlike in Fourier space, the fluctuations behave very differently in two and three dimensions). Fur- 
thermore, in these simulations we do not include a stochastic flux in the concentration equation, 
i.e., we set * = 0, so that all fluctuations in the concentration arise from being out of thermody- 
namic equilibrium. With this approximation we do not need to model the chemical potential of 
the mixture and obtain fj.^. It is known experimentally that the nonequilibrium fluctuations are 
much larger than the equilibrium ones for the conditions we consider [30]. 



In the simple linearized theory presented in Section II D 2 several approximations are made. 
The ffist one is that a quasi-periodic approximation is used even though the actual system is not 
periodic along the y direction. This source of error has already been studied numerically in Ref. 
|22j . A second source of error is the use of a Boussinesq approximation. In this approximation, it 
is assumed that pi = p + Ap/2 and P2 = p — A/9/2, where Ap is a small density difference between 
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the two fluids, Ap/p <c 1, so that density is approximately constant and /3 <c 1. More precisely, in 
the Boussinesq model the gravity term in the velocity equation only enters through the product 
Pg so the approximation consists of taking the limit (3^0 and g ^ oo while keeping the product 
f3g fixed. The final approximation made in the simple theory is that the transport coefficients, i.e., 
the viscosity and diffusion coefficients, are assumed to be constant. Here we evaluate the validity 
of the constant-coefficient constant-density approximation (p, rj and x constant, /3 ^ O), as well as 
the constant-density (Boussinesq) approximation alone {p constant, /3 -> 0, but variable 77, by 
comparing with the solution to the complete low Mach number equations (p, 77, x P variable). 



A. Simulation Parameters 



We base our parameters on the experimental studies of diffusive mixing in a water-glycerol 
mixture, as reported in Ref. [40]. The physical domain is 1cm x 0.25 cm discretized on a uniform 
128 X 32 two dimensional grid, with a thickness of 1 cm along the z direction. Gravity is applied 



along the negative y (vertical) direction. Reservoir boundary conditions ( 13 ) are applied in the y- 
direction and periodic boundary conditions in the a;-direction. We set the concentration to c = 0.39 
on the bottom boundary and c = on the top boundary, and apply no-slip boundary conditions 
for the velocity at both boundaries. The initial condition is c(t = 0) = 0.39(y/0.25 — 1), which is 
close to the deterministic steady-state profile. A very good fit to the experimental equation of 
state (dependence of density on concentration) over the whole range of concentrations of interest is 



provided by the EOS (11) with the density of water set to p2 = Ig/cm^ and the density of glycerol 
set to pi — 1.29 g/cm^. In these simulations the magnitude of the velocity fluctuations is very small 
and we did not use filtering (see Appendix [A|) . 

Experimentally, the dependence of viscosity on glycerol mass fraction has been fitted to an 
exponential function [40], which we approximate with a quadratic function over the range of con- 
centrations of interest, 

v{c) = exp(2.06c + 2.32c'') w i/o(1.0 + 0.66c + 12c'), (52) 
where experimental measurements estimate i^o ~ lO^'cm'/s. The diffusion coefficient dependence 
on the concentration is, to our knowledge, not known, and is in fact strongly affected by thermal 
fluctuations and spatial confinement (TJ [50] . We approximate the dependence assuming a Stokes- 
Einstein relation, 

x(c) = ^«Xo(1.0-2.2c+1.2c'), (53) 
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where experimental estimates for water-glycerol mixtures give xo ~ 10"^cm^/s, with a Schmidt 
number Sc — J^/x~ 10'^- This very large separation of scales between mass and momentum diffusion 
is not feasible to simulate with our explicit temporal integration methods. Referring back to the 



simplified theory (25), we see that the shape of the spectrum of the steady-state concentration 
fluctuations, and in particular, the cutoff wavenumber due to gravity, is determined from the 
product and not x and individually. Therefore, as also done in Ref. ^22], we choose xo 
and 1^0 so that x (c) (c) is kept at the physical value of 10^'' (cm^/s)^ but the Schmidt number is 
reduced by two orders of magnitude, — iy{c) /x(c) = 10, where c = 0.39/2 is an estimate of the 
average concentration. The condition v (c) ~ 10~^ cm^/s and x (c) ~ 10"" cm^/s gives our simulation 
parameters i^o « 6.1 x 10"*cm^/s and xo ~ 1-6 x lO^^cm^/s. 

The physical value for gravity is 5 ~ 10^ cm/s^ and the solutal expansion coefficient (3 (c) « 0.234 
follows from pi and p2. When employing the Boussinesq approximation, in which gravity only 
enters through the product /3g, we set pi — 1.054 and p2 = 1.044 so that (3 — 0.01 and increase 
gravity by the corresponding factor to g = 2.34 • 10" cm/s^ in order to keep fixed at the physical 
value. We also performed simulations with a weaker gravity, g « lO^cm/s^, which enhances the 
giant fluctuations. 

B. Results 

We employ the explicit midpoint temporal integrator (which we recall is third-order accurate for 
static covariances) and set At — 0.005 s, which results in a diffusive Courant number vlS.t/ /S.x'^ « 0.1. 
We skip the first 50,000 time steps (about 5 diffusion crossing times) and then collect samples from 
the subsequent 50,000 time steps. We repeat this eight times to increase the statistical accuracy 



and estimate error bars. To compare to the theory (25), we set the concentration gradient to /ly = 
0.39/0.25 cm~^ and evaluate p ~ 1.05^ at c = 0.39/2 from the equation of state. When computing 
the theory, we account for errors in the discrete approximation to the continuum Laplacian by 
using the effective wavenumber 

sin (fc,Aa;/2) 

^^-^^ (fc.Ax/2) ^^^> 

instead of the actual discrete wavenumber [22] . 

The results for the static spectrum of concentration fluctuations S^^^ {k^,ky = 0) = (('^'^) {^'^ 
as a function of the modified wavenumber fc^ (54) are shown in Fig. [2j When there is no gravity, we 



see the characteristic giant fluctuation power-law spectrum of the fluctuations, modulated at small 
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Figure 2: Comparison between the simple theory (25) (hnes) and numerical results (symbols) in the ab- 
sence (black circles) of gravity and with standard gravity g ~ 10^ ^ present (the cutoff wavenumber 
kg K, 246 cm^^), for the complete variable-coefficient variable-density low Mach model (green triangles), 
the variable-coefficient constant-density (Boussinesq) approximation (blue diamonds), and the constant- 
coefficient constant-density approximation (red squares). Also shown are results for a weaker gravity, 
g sa 10^ ^ (the cutoff wavenumber kg f38 cm^^), for the complete variable-coefficient variable-density 
low Mach model (magenta pluses) and the constant-coefficient constant-density approximation (cyan stars). 



wavenumbers due to the presence of the physical boundaries [22j. When gravity is present, fluc- 
tuations at wavenumber below the cutoff kg — [h^^gji/ {vx)]''''^ are suppressed. If we use a constant- 
coefficient approximation, in which we reduce /? — O.Of so that p p{c) and also fix the transport 
coefficients at v{c) = v{c) and x(c) = x(c), we observe good agreement with the quasi-periodic 



theory (25). However, if we make the transport coefficients dependent on the concentration as 



m ( [52|53l ), we see a substantial deviation from the simplified theory. If we further remove the 
Boussinesq approximation and employ the physical dependence of density on concentration, we 
see little change. This shows that under the sort of parameters present in the experiments on 
diffusive mixing in water-glycerol mixture, it is reasonable to make the Boussinesq incompressible 
approximation; however, the spatial dependence of the viscosity and diffusion coefficient cannot be 
ignored if quantitative agreement is desired. Even though the constant-coefficient approximation 
gives qualitatively the correct shape and a better choice of the constant transport coefficients may 
improve its accuracy, there is no obvious or simple procedure to a priori estimate what parameters 
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should be used. 

VI. DIFFUSIVE MIXING IN HARD-DISK AND HARD-SPHERE FLUIDS 

In this section we study the appearance of giant fluctuations during time- dependent diffusive 
mixing. As a validation of the low Mach number fluctuating equations and our algorithm, we 
perform simulations of diffusive mixing of two fluids of different densities in two dimensions. We 
find excellent agreement between the results of low Mach number (continuum) simulations and 
hard-disk molecular dynamics (particle) simulations. This nontrivial test clearly demonstrates the 
usefulness of low Mach number models as a coarse-grained mesoscopic model for problems where 
sound waves can be neglected. 

Our simulation setup is illustrated in Fig. [3} We consider a periodic square box of length L 
along both the x (horizontal) and y (vertical) directions, and initially place all of the fluid of species 
one (colored red) in the middle third of the domain, i.e., we set c — I for L/3 < y < 2L/3, and 
c = otherwise, as shown in the top left panel of the figure. The two fluids mix diffusively and 
at the end of the simulation the concentration field shows a rough diffusive interface as confirmed 
by molecular dynamics simulations shown in the top right panel of the figure. The deterministic 
equations of diffusive mixing reduce to a one dimensional model due to the translational symmetry 
along the x axes, and would yield a flat diffusive interface as illustrated in the bottom left panel 
of the figure. However, fluctuating hydrodynamics correctly reproduces the interface roughness, as 
illustrated in the bottom right panel of the figure and demonstrated quantitatively below. 

A. Molecular Dynamics 

We performed Hard Disk Molecular Dynamics (HDMD) simulations of diffusive mixing using 
a modification of the public-domain code developed by the authors of Ref. |59j . We use arbitrary 
(molecular) units of length, time and mass for convenience. All hard disks had a diameter a — 1, 
and we set the temperature at fc^T = 1. The molecular mass for the first fluid component was fixed 
at TOi = 1, and for the second component at m2 — Rrrii. For mass ratio R — 1, the two types of 
disks are mechanically-identical and therefore the species label is simply a red-blue coloring of the 
particles. In this case p2 = pi and the low Mach number equations reduce to the incompressible 
equations of fluctuating hydrodynamics with a passively-advected concentration field. For the case 
of unequal particle masses, mechanical equilibrium is obtained if the pressures in the two fluid 
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Figure 3: Diffusive mixing between two fluids of unequal densities, R — P2lpi—^, with coloring based on 
concentration, red for the pure first component, c — \, and blue for the pure second component, c = 0. 
A smoothed shading is used for the coloring to eliminate visual discretization artifacts. The simulation 
domain is periodic and contains 128^ hydrodynamic (finite volume) cells. The top left panel shows the 
initial configuration, which is the same for all simulations reported here. The top right panel shows the 
final configuration at time t = 5, 800 as obtained using molecular dynamics. The bottom left panel shows 
the final configuration obtained using deterministic hydrodynamics, while the right panel shows the final 
configuration obtained using fluctuating hydrodynamics. 

components are the same. It is well-known from statistical mechanics that for hard disks the 
pressure is 

P = y ((/))• n • kgT 

where n = N/V is the number density, and Y (0) is a prefactor that only depends on the packing 
fraction (fraction of the plane occupied by the disks) <j) — n (na'^/i) and not on the molecular mass. 
Therefore, for a mixture of disks with equal diameters, at constant pressure, the number density 
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and the packing fraction tp are constant independent of the composition. The equation of state at 
constant pressure and temperature is therefore 

n n rnrii nm2 ' 



which is exactly of the form (11) with pi = nmi and p2 — rim2. The chemical potential of such a 



mixture has the same concentration dependence as for a low-density gas mixture 

P-^^ksT = c (1 — c) [cm-i + (1 — c) 7712] . 

We used a packing fraction of = 0.6 for all simulations reported here. This packing fraction 
is close to the freezing transition point but is known to be safely in the (dense) gas phase (there 
is no liquid phase for a hard-disk fluid). The initial particle positions were generated using a 
nonequilibrium molecular dynamics simulation as in the hard-particle packing algorithm described 
in Ref . [60] . After the initial configuration was generated the disks were assigned a species according 
to their y coordinate, and the mixing simulation performed using event-driven molecular dynamics. 

In order to convert the particle data to hydrodynamic data comparable to that generated by 
the fluctuating hydrodynamics simulations, we employed a grid of hydrodynamic cells that were 
each a square of linear dimension = 10. At the chosen packing fraction — 0.6 this corresponds 
to about 76 disks per hydrodynamic cell, which is deemed a reasonable level of coarse-graining for 
the equations of fluctuating hydrodynamics to be a reasonably-accurate model, while still keeping 
the computational demands of the simulations manageable. We performed HDMD simulations for 
systems of size — 64 and = 128 cells, and simulated the mixing process to a final simulation 
time of t = 5, 800. The largest system simulated had about 1.25 million disks (each simulation took 
about 5 days of CPU time), which is well into the "hydrodynamic" rather than "molecular" scale. 

Every 58 units of time, particle data was converted to hydrodynamic data for the purposes 
of analysis and comparison to hydrodynamic calculations. There is not a unique way of coarse- 
graining particle data to hydrodynamic data |6H 162] : however, we believe that the large-scale 
(giant) concentration fluctuations we study are not affected by the particular choice. We therefore 
used a simple method consistent with the philosophy of finite- volume conservative discretizations. 
Specifically, we coarse-grained the particle information by sorting the particles into hydrodynamic 
cells based on the position of their centroid, as if they were point particles. We then calculated pi 
and P2 in each cell based on the total mass of each species contained inside the given cell. Since all 
particles have equal diameter other definitions that take into account the particle shape and size 
give similar results. 
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1. Viscosity v 

The hydrodynamic simulations described below require as input transport coefficients, notably, 
the shear viscosity 77 and diffusion coefficient x- As discussed in more detail in Ref. [SD], these coef- 
ficients are not universal material constants but rather depend on the spatial scale under question. 
We emphasize that this scale-dependent renormalization is not a molecular scale effect but rather 
an effect arising out of hydrodynamic fluctuations, and persists even at the hydrodynamic scales we 
are examining here. The best way to define and measure transport coefficients is by examining the 
dynamics of equilibrium fluctuations, specifically, by examining the dynamic structure factors of 
the hydrodynamic fields [I6j, i.e., the equilibrium averages of the spatio-temporal Fourier spectra of 
the fluctuating hydrodynamic fields. For a hydrodynamic variable x that is transported by a purely 
diffusive process, the spectrum of the fluctuations at a given wavenumber k and wavefrequency oj 
is expected to be a Lorentzian peak of the form 

S,{k,uj) = {x{k,u)x* {k,uj)) - [u' + CkX\ 
where in general the diffusion constant C (fc) depends on the the wavenumber k (wavelength A = 
27r/fc). We can therefore estimate the diffusion coefficient x by fltting a Lorentzian peak to {k,Lo) 
for different fc's (i.e., x = c). Similarly, we can estimate the kinematic viscosity u = 77/p by fltting a 
Lorentzian curve to dynamic structure factors for the scaled vorticity, x = k^^ (V x v)_,. 

We performed long equilibrium molecular dynamics simulations of systems corresponding to 
a grid of — 32 hydrodynamic cells, and then calculated the discrete spatio-temporal Fourier 
spectrum of the hydrodynamic flelds at a collection of discrete wavenumbers k. Since these simu- 
lations are at equilibrium, the systems are well-mixed, speciflcally, the initial conflgurations were 
generated by randomly assigning a species label to each particle. We then performed a nonlinear 
least squares Lorentzian flt in lo for each fc and estimated the width of the Lorentzian peak. The 
results for the dynamics of the equilibrium vorticity fluctuations are shown in Fig. |4j We see 
that kinematic viscosity is relatively constant for a broad range of wavelengths, consistent with 
fluctuating hydrodynamics calculations [63] and previous molecular dynamics simulations [64J. For 
the pure component one fluid, c = 1, with density p ~ 0.764 the figure shows v « 3.3. We therefore 
used ?7i w 0.764 • 3.3 ~ 2.5 in all of the hydrodynamic runs reported below. This is about 20% higher 
than the prediction of the simple Enskog kinetic theory [65], r] « 2.06, and is consistent with the 
estimates reported in Ref. |64j . Because of the smallness of the diffusion coefficient at the densities 
we study, more specifically, because the Schmidt number = v/x is larger than 10, we were unable 
to obtain reliable estimates for x {k) from the dynamic structure factor for concentration. 
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Figure 4: Estimates of the momentum diffusion coefficient (viscosity) v — rj/ p obtained from the width of 
the central peak in the dynamic structure factor of vorticity. A collection of 24 distinct discrete wavenumbers 
fe were used and the width of the peaks estimated using a nonlinear least squares Lorentzian fit. 



Simple dimensional analysis or kinetic theory shows that -q 
the pure second fluid component is 



% = ■ni\i — = vi^- 



m, and therefore the viscosity of 



(55) 



There is no simple theory that accurately predicts the concentration dependence of the viscosity of 
a hard disk mixture at higher densities [66j . To our knowledge there is no published Enskog kinetic 
theory calculations for hard-disk mixtures in two dimensions, even for the simpler case of equal 
diameters. As an approximation to the true dependence, we employed a simple linear interpolation 
of the kinematic viscosity v{c) — 7]{c) /p as a function of the mass concentration c between the two 
known values = v {c ~ 1) ^ 3.3 and — {c = 0) ^ v^j^fR. The numerical results for mixtures 
with mass ratios R = 2 and i? = 4 in Fig. |4] are consistent with this approximation to within the 
large error bars. For example, for c— 1/2 and R = 4 the interpolation gives v = 3 - 3.3/4 w 2.5 which 
is in reasonable agreement with the numerical estimate. 
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B. Hydrodynamics 

We now turn to hydro dynamic simulations of the diffusive mixing of hard disks. Our hydrody- 
namic calculations use the same grid of cells used to convert particle to hydrodynamic data. The 
only input required for the hydrodynamic calculations, in addition to those provided by equilib- 
rium statistical mechanics, are the transport coefficients of the fluid as a function of concentration, 
specifically, the shear viscosity 77, the diffusion coefficient x, and the bulk viscosity. In our hy- 
drodynamic calculations we ignore the bulk viscosity term in the shear stress. Given that in the 
linearized theory the bulk viscosity drops out and that it is quite difficult to measure in particle 
methods or in experiments, this approximation is expected to not affect the results significantly, 
as will be justified a posteriori by the good match. 

The values for the transport coefficients used in the spatio-temporal discretization, as explained 
in Ref. [50], are not material constants independent of the discretization. Rather, they are bare 
transport values ryo and xo measured at the length scales of the grid size. We assumed that the 
bare transport coefficients obey the same scaling with the mass ratio R as predicted by Enskog 
kinetic theory ( 55|56 ). As explained earlier, theoretical arguments and molecular dynamics results 



suggest that renormalization effects for viscosity are small and can be safely neglected. We have 
therefore fixed the viscosity in the hydrodynamic calculations based on the molecular dynamics 



estimate rjo — 2.5 for the pure fiuid with molecular mass m = 1 (see Section VIAl). However, the 
bare diffusion coefficient is strongly dependent on the size of the hydrodynamic cells (held fixed in 
our calculations at Ax = Ay = 10), and on whether filtering (see Section |A]) is used. Therefore, the 
value of Xo needs to be adjusted based on the spatial discretization, in such a way as to match the 
behavior of the molecular dynamics simulations at length scales much larger than the grid spacing. 
We describe the exact procedure we used to accomplish this in the next section. 

The time step in our explicit algorithm is limited by the viscous CFL number — vAt/Ax^ < 
1/4. Since the hydrodynamic calculations are much faster compared to the particle simulations, 
we used the more expensive RK3 temporal integrator with a relatively small time step At — 1.45, 
corresponding to ~ 0.05 for c — 1. For R ~ 1 and — 64 we employed a larger time step. 
At = 3.625 {a^ w 0.125), with no measurable temporal discretization artifacts for the quantities we 
study here. We are therefore confident that the discretization errors in this study are dominated 
by spatial discretization artifacts. In future work we will explore semi-implicit discretizations and 
study the effect of taking larger time steps on the temporal accuracy. Note that at these parameters 
for c = 1 the isothermal speed of sound is ct ~ 5.1 so that a compressible scheme would require a 
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time step on the order of Ai ^ 1 (corresponding to advective CFL of about a half). By contrast, 
the exphcit low Mach number algorithm is stable for At < 7.5. This modest gain is due to the 
small hydrodynamic cell we use here in order to compare to molecular dynamics. For mesoscopic 
hydrodynamic cells the gain in time step size afforded by the low Mach formulation will be several 
orders of magnitude larger. 

For mass ratio R — 1 and R = 2, the hydrodynamic calculations were initialized using statistically 
identical configurations as would be obtained by coarse-graining the initial particle configuration. 
This implies a sharp, step-like jump in concentration at y = L/3 and y ~ 2L/?,. Since our spatio- 
temporal discretization is not strictly monotonicity-preserving, such sharp concentration gradients 
combined with a small diffusion coefficient xo lead to a large cell Peclet number. This may in turn 
lead to large deviations of concentration outside of the allowed interval < c < 1 for larger mass 
ratios. Therefore, for i? = 4 we smoothed the initial condition slightly so that the sharp jump in 
concentration is spread over a few cells, and also employed a 9 point filter {wp = 4, see Appendix 
[A| ). We verified that for R — 2 using filtering only affects the large wavenumbers and does not 
appear to affect the small wavenumbers we study here, provided the bare diffusion coefficient xo is 
adjusted based on the specific filtering width wp- 



1. Diffusion Coefficient x 

For the inter-species diffusion coefficient x, which we emphasize is distinct from the self-diffusion 
coefficients for particles of either species, Enskog kinetic theory predicts no concentration depen- 
dence and a simple scaling with the mass ratio |66j , 

X(i?) = x(i? = l)/^. (56) 
This particular dependence on mass ratio R comes from the fact that the average relative speed 
between particles of different species is ^ ^fcaT/m^, where m^j — mimi/ (rrii + ma) is the reduced 
molecular mass. We have assumed in our hydrodynamic calculations that the diffusion coefficient 



is independent of concentration and follows ( 56 ) . The only input to the hydrodynamic calculation 
is the bare self-diffusion coefficient for the pure component fluid, xo (R = !)• Diffusion is strongly 
renormalized by thermal fluctuations, and fluctuating hydrodynamics theory and simulations pre- 
dict a strong dependence of the diffusion coefficient x on the wavelength \50\, consistent with 
molecular dynamics results 

In order to estimate the appropriate value of the bare diffusion coefficient xo we numerically 
solved an inverse problem. Using simple bisection, we looked for the value of xo that leads to best 
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Figure 5: (Left panel) Diffusive evolution of the horizontally-averaged density p^''-' (y) for a system of size 
Nc — 64 hydrodynamic cells and density ratio i? = 1, as obtained from HDMD simulations (circles, averaged 
over 64 runs), deterministic hydrodynamics with XcS = 0.2 (dashed lines), and fluctuating hydrodynamics 
with xo = 0.09 (squares, averaged over 64 runs). Error bars are comparable to the symbol size and not 
shown for clarity. (Right panel) Same as the left panel except the density ratio is i? = 2 and the transport 
coefficients are adjusted according to ( 55|56 1. 



agreement for the average or "macroscopic" diffusion (mixing) between the particle and continuum 
simulations. Specifically, we calculated the density of the first species pi''' (y) along the y-direction 
by averaging pi in each horizontal row of hydrodynamic cells. The results for pi'"' for mass ratios 
R — 1 and R — 4 are shown in Fig. [5] at different points in time for systems of size = 64 cells. 
The figures show the expected sort of diffusive mixing profile, and is exactly what would be used 
in experiments to measure diffusion coefficients using fluorescent techniques such as Fluorescence 
Recovery After Photo-bleaching (FRAP) [67J. This macroscopic measurement smooths over the 
fluctuations (roughness) of the diffusive interface and only measures an effective diffusion coefficient 
at the scale of the domain length L. If deterministic hydrodynamics is employed, pi'"' (y) is the 
solution of a one-dimensional system of equations obtained by simply deleting the stochastic forcing 
and the x-dependence in the low Mach equations. Instead of solving this system analytically, we 
employed our spatio-temporal discretization with fiuctuations turned off, and with an effective 
diffusion coefficient x = Xcs that accounts for the renormalization of the diffusion coefficient by the 
thermal fluctuations. 

By matching the profile pi'"' (y) between the HDMD and the fluctuating and deterministic hy- 
drodynamic simulations at mass ratio R — 1 and system size — 64 cells, we obtained estimates 
for the bare xo and the renormalized Xos coefficient (see Fig. [5|. The best estimate for the bare 
diffusion coefficient based on this matching in the absence of filtering is xo = 0.09 ± O.Of . This 
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compares reasonably-well to the prediction of Enskog theory [65] of x ~ 0.08, as well as to the 

measurement of the self-diffusion coefficient for a periodic system with 169 disks reported in Ref. 
j64| . X ~ 0.14 (recall that a single hydrodynamic cell in our case contains about 76 particles). When 
a 5-point filter is employed the estimate is xo {wf ~ 2) ^ 0.12 and when a 9-point filter is employed 
Xo {wp = 4) w 0.11. The estimated renormalized diffusion coefficient is much larger, Xcs ~ 0.20±0.01, 
consistent with a rough estimate based on the simple theory presented in Ref. \50{ . 



To within statistical accuracy we were not able to detect the increase in the estimated diffusion 
coefficients when using the larger systems of size — 128 cells, however, for — 32 it was clear 
that Xcs is reduced. In any case it should be noted that the theory in Ref. [50] is for steady-state 
diffusion and does not strictly apply in the present context. 

It is important to emphasize that Xcs is not a material constant but rather depends on the details 
of the problem in question, in particular, the system geometry and size, the initial and boundary 
conditions, and probably also the time at which the measurement is performed. By contrast, xo is a 
constant for a given spatial discretization, and one can use the same number for different scenarios so 
long as the hydrodynamic cell size and the filter are kept fixed. Unlike deterministic hydrodynamics, 
which presents an incomplete picture of diffusion, fiuctuating hydrodynamics correctly accounts for 
the important contribution of the thermal velocity fiuctuations and the roughness of the diffusive 
interface seen in Fig. [3] 



In order to compare the molecular dynamics and the hydrodynamic simulations we calculated 
several statistical quantities: 

1. The averages of pi along the directions perpendicular to the concentration gradient. 



where the integral is discretized as a direct sum over the hydrodynamic cells. Note that it 
is statistically better to use conserved quantities for such macroscopic averages than to use 
non-conserved variables such as concentration [68]. 

2. The concentration averaged along the direction of the gradient, 




for = 64 



for TV, = 128 



C. Comparison between particle and continuum simulations 





y=0 
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Intuitively, this is a measure of the thickness of the red strip in Fig. [3j 

3. The y-coordinate of the "center-of-mass" of concentration along the direction perpendicular 
to the gradient, 

{x) = L'^ / y-c{x,y) dy. 
Jy=o 

Intuitively, this is a measure of the height of the centerline of the red strip in Fig. [Sj 

All quantities were sampled at certain pre-specified time points in a number of statistically- 
independent simulations and then means and standard deviations calculated from the data 
points. For system of size = 64 cells we used — 64 simulations, and for systems of size ~ 128 
we used ~ 32 simulations. By far the majority of the computational cost was in performing the 
HDMD simulations. 



1. Average Concentration Profiles 

Once xo and Xcff were estimated based on simulations of a constant-density {R = l) fluid (see 
Section [VI B 1 ) , scaling theory ( 55|56 ) can be used to estimate them for different density ratios. 



In Fig. [5] we show p[^^ (y) for mass ratio R — 2, showing good agreement between HDMD and 
hydrodynamics, especially when fluctuations are accounted for. For R = i a direct comparison is 
difficult because the initial condition was slightly different in the hydrodynamic simulations due to 
the need to smooth the sharp concentration gradient for numerical reasons, as explained earlier. 
This difference strongly affects the shape of pj''' (y) at early times, however, it does not significantly 
modify the roughness of the interface, which we study next. 



2. Interface Roughness 

The most interesting contribution of fluctuations to the diffusive mixing process is the ap- 
pearance of giant concentration fluctuations in the presence of large concentration gradients, as 
evidenced in the roughness of the interface between the two fluids during the early stages of the 
mixing. In order to quantify this interface roughness we used the one-dimensional power spectra 

Sc {K) = {cvcD and S,, {K) = (hjil) . 
Note that here we do not correct the discrete wavenumber for the spatial discretization artifacts 
and continue to use instead of k^. In a time-dependent context the discretization errors are not 
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Figure 6: Discrete spatial spectrum of the interface fluctuations for R = 1 and Nc = 128 at several points 
in time (averaged over 32 simulations), for fluctuating hydrodynamics (squares with error bars) and HDMD 
(circles, error bars comparable to those for squares). Note that the largest wavenumber supported by the grid 
is fcmax = tt/Ax ~ 0.314. The larger wavenumbers are however dominated by spatial truncation errors and 
the filter employed (if any) and we do not show them here. (Left panel) Spectrum Sc (k^) of the vertically- 
averaged concentration. (Right panel) Spectrum Sh (kx) of the position of the vertical "center-of-mass" of 
concentration. 
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Figure 7: Same as Fig. |6]but for density ratio R — A. 



known in general, and, in any case, our main goal is to compare particle and continuum simulations, 
rather than compare to theory. 

The temporal evolution of the spectra 5c and Su is shown in Fig. [6] for mass ratio i? = 1, and 
in Fig. [7] for mass ratio i? = 4, for both HDMD and low Mach number fluctuating hydrodynamics 
(note that deterministic hydrodynamics would give identically zero for any spectral quantity). 
We observe an excellent agreement between the two, including the correct initial evolution of the 
interface fluctuations. 
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Figure 8: Mixing to a time 128 times longer than previous results, with results reported at time intervals 
t — 7424 for i = 1, . . . , 10. These long simulations are only feasible for the fluctuating hydrodynamics 
code, and employ a somewhat larger time step At = 3.625. {Left) Horizontally-averaged pi, as shown for the 
shorter runs in the left panel of Fig. [s] (Right) The spectnun of interface fluctuations 5*^ (k^), as shown in 
the left panels of Figs. |6]and[7]for the shorter runs. The theoretical estimates for the spectrum of equilibrium 
fluctuations, which is independent of wavenumber, is also shown. We also indicate the theoretical prediction 
for the power-law of the spectrum of steady-state nonequilibrium fluctuations under an applied concentration 
gradient, Sc ^ fc~*. 

Note that for a finite system, eventually complete mixing will take place and the concentration 
fluctuations will have to revert to their equilibrium spectrum, which is flat in Fourier space instead 
of the power-law behavior seen out of equilibrium. In Fig. [8] we show results for mixing up to a 
time t — 7.42 • 10^ (this is 128 times longer than those described above). These long simulations 
are only feasible for the fluctuating hydrodynamics code, and employ a somewhat larger time 
step At — 3.625. The results clearly show that at late times the spectrum of the fluctuations 
reverts to the equilibrium one; however, this takes some time even after the mixing is essentially 
complete. Linearized incompressible fluctuating hydrodynamics \22\ [36] predicts that at steady 
state the spectrum of nonequilibrium concentration fluctuations is a power law with exponent —4, 

^ (Vc)^ fc"*. The dynamically-evolving spectra in the right panel of Fig. [sjshow approximately 
such power-law behavior for intermediate times and wavenumbers. 

D. Three Dimensions 

In order to illustrate the appearance of giant fluctuations in three dimensions we performed 
simulations of mixing in a mixture of hard spheres with equal diameters, a — I, and mass ratio 
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R~ A. The packing density was chosen iohe cf) — 0.45, which corresponds to a very dense gas, but is 
still well below the freezing point = 0.49. For the hydrodynamic simulations we used cubic cells 
of dimension Ax = 5, which corresponds to about 107 particles per hydrodynamic cell on average. 
In Fig. [9] we show results from a single simulation with a grid of size 128 x 64 x 128 cells, which 
would correspond to about 10* particles. This makes molecular dynamics simulations infeasible, 
and makes hydrodynamic calculations an invaluable tool in studying the mixing process at these 
mesoscopic scales. 

In the hydrodynamic simulations we used bare transport coefficient values based on Enskog 
kinetic theory for the hard-sphere fluid [HS]- For the single-component fluid with molecular mass 
m = 1, this theory gives % ~ 2.32 and xo ~ 0.053, which corresponds to a bare Schmidt number 
Sc — t'o/xo ~ 51. We employed the same model dependence of bare transport coefficients on 
concentration as for hard disks, see Eqs. ( 55|56 ). The time step was set at At = 1 (corresponding 
to viscous CFL number /? = I'oAt/Ax'^ ~ 0.1). In three dimensions, the cell Peclet number is reduced 
with decreasing Ax and we did not find it necessary to employ any filtering. 

Instead of the fully periodic domain used in the two dimensional hard-disk simulations, here we 



employ the fixed-concentration boundary conditions ( 13 ) and set c{y = 0; t) = at the bottom and 
c{y — Ly] t) — 1 at the top boundary. This emulates the sort of "open" or "reservoir" boundaries 
|69j that mimic conditions in experimental studies of diffusive mixing |40j . The initial condition 
is a fully phase-separated mixture with c = 1 for y > L/2, and c ~ otherwise. As the mixing 
process continues the diffusive interface roughens and giant concentrations appear, as illustrated 
in Fig. [9] and also observed experimentally in water-glycerol mixtures in Ref. |40j . In three 
dimensions, however, the diffusive interface roughness is much smaller than in two dimensions, 
being on the order of only 20 molecular diameters for the snapshot shown in the figure. This 
illustrates the importance of dimensionality when including thermal fluctuations. In particular, 
unlike in deterministic fluid dynamics, in fluctuating hydrodynamics one cannot simply eliminate 
dimensions from consideration even in simple geometries. 

Approximate theory based on the Boussinesq approximation and linearization of the equations 
of fluctuating hydrodynamics has been developed in Ref. [36j and applied in the analysis of ex- 
perimental results on mixing in a water-glycerol mixture in the presence of gravity [ID]. The 
simulations reported here do not make the sort of approximations necessary in analytical theories 
and can in principle be used to study the mixing process quantitatively. However, it is important 
to emphasize that in realistic liquids, such as a water-glycerol mixture, the Schmidt number is on 
the order of a thousand. This makes explicit time stepping schemes that fully resolve the dynamics 
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Figure 9: Diffusive mixing in three dimensions similar to that illustrated in Fig. [3] for two dimensions. 
Parameters are based on Enskog kinetic theory for a hard-sphere fluid at packing fraction (j) — 0.45, and 
there is no gravity. The mixing starts with the top half being one species and the bottom half another 
species, with density ratio R = A, and concentration is kept fixed at the top and bottom boundaries while 
the side boundaries are periodic. A snapshot taken at time t = 5, 000 is shown. {Top panel) The side panes 
show two dimensional slices for the concentration c. The approximated contour surface c = 0.2 is shown with 
color based on surface height to illustrate the rough diffusive interface. (Bottom left panel) Similar as top 
panel but bottom pane shows vertically-averaged concentration (x, z), illustrating the giant concentration 
fluctuations. (Bottom right panel) The Fourier spectrum Sc {kx, ky) of c„. The color axes is logarithmic and 
clearly shows the appearance of large scale (small wavenumber) fluctuations, as also seen in Fig. [T] in two 
dimensions. 

of the velocity fluctuations infeasible. In future work we will consider semi-implicit type stepping 
methods that relax the severe time stepping restrictions present in the explicit schemes considered 
here. 
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VII. CONCLUSIONS 

The behavior of fluids is strongly affected by thermal fluctuations at scales from the micro- 
scopic to the macroscopic. Fluctuating hydrodynamics is a powerful coarse-grained model for fluid 
dynamics at mesoscopic and macroscopic scales, at both a theoretical and a computational level. 
Theoretical calculations are rather complicated in the presence of realistic spatial inhomogeneities 
and nontrivial boundary conditions. In numerical simulations, those effects can readily be handled, 
however, the large separation of time scales between different physical processes poses a fundamen- 
tal difficulty. Compressible fluctuating hydrodynamics bridges the gap between molecular and 
hydrodynamic scales. At spatial scales not much larger than molecular, sound and momentum and 
heat diffusion occur at comparable time scales in both gases and liquids. At mesoscopic and larger 
length scales, fast pressure fluctuations due to thermally-actuated sound waves are much faster 
than diffusive processes. It is therefore necessary to eliminate sound modes from the compressible 
equations. In the deterministic context this is accomplished using low Mach number asymptotic 
expansion. 

For homogeneous simple fluids or mixtures of dynamically-identical fluids the zcroth order low 
Mach equations are the well-known incompressible Navier-Stokes equations, in which pressure is a 
Lagrange multiplier enforcing a divergence-free velocity field. In mixtures of dissimilar fluids, local 
changes in composition and temperature cause local expansion and contraction of the fluid and thus 
a nonzero velocity divergence. In this paper we proposed low Mach number fluctuating equations 
for isothermal binary mixtures of incompressible fluids with different density, or a mixture of low- 
density gases with different molecular masses. These equations are a straightforward generalization 
of the widely-used incompressible fluctuating Navier-Stokes equations. In the low Mach number 
equations the incompressibility constraint V • v = is replaced hy V ■ v = —l3{Dc/Dt), which 
ensures that compositional changes are accompanied by density changes in agreement with the fluid 
equation of state (EOS) at constant pressure and temperature. This seemingly simple generalization 
poses many non-trivial analytical and numerical challenges, some of which we addressed in this 
paper. 

At the analytical level the low Mach number fluctuating equations are different from the incom- 
pressible equations because the velocity divergence is directly coupled to the time derivative of the 

concentration fluctuations. This means that at thermodynamic equilibrium the velocity is not only 
white in space, a well-known difficulty with the standard equations of fluctuating hydrodynamics, 
but is also white in time, adding a novel type of difficulty that has not heretofore been recognized. 
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The unphysically fast fluctuations in velocity are caused by the unphysical assumption of infinite 
separation of time scales between the sound and the difi^usive modes. This unphysical assumption 
also underlies the incompressible fluctuating Navier-Stokes equations, however, in the incompress- 
ible limit /? — 7> the problem is not apparent because the component of velocity that is white in 
time disappears. Here we analyzed the low Mach equations at the linearized level, and showed that 
they reproduce the slow diffusive fluctuations in the full compressible equations, while eliminating 
the fast pressure fluctuations. At the formal level, we suggest that a generalized Hodge decom- 
position can be used to separate the vortical (solenoidal) modes of velocity as the independently 
fluctuating variable, and a gauge formulation used to eliminate the divergence constraint. Such 
nonlinear analysis is deferred for future research, and here we relied on the fact that the temporal 
discretization regularizes the short-time dynamics at time scales faster than the time step size At. 

At the numerical level, the low Mach number equations pose several distinct challenges. The 
flrst challenge is to construct conservative spatial discretizations in which density is advected in 
a locally-conservative manner while still maintaining the equation of state constraint relating the 
local densities and composition. We accomplish this here by using a specially-chosen model EOS 
that is linear yet still rather versatile in practice, and by advecting densities using a velocity that 
obeys a discrete divergence constraint. As in incompressible hydrodynamics, enforcing this con- 
straint requires a nonlocal Poisson pressure solver that dominates the computational cost of the 
algorithm. A second challenge is to construct temporal integrators that are at least second-order 
in time. We accomplish this here by formally introducing an unconstrained gauge formulation of 
the equations, while at the same time taking advantage of the gauge degree of freedom to avoid 
ever explicitly dealing with the unphysical gauge variable. The present temporal discretizations 
are purely explicit and are similar in spirit to an explicit projection method. A third and re- 
maining challenge is to design efficient temporal integrators that handle momentum diffusion, the 
second- fastest physical process, semi-implicitly. This poses well-known challenges even in the in- 
compressible setting. These challenges were bypassed in recently-developed temporal integrators 
for the incompressible fiuctuating Navier-Stokes equations [22] by avoiding the splitting inherent 
in projection methods. Extending this type of Stokes-system approach to the low Mach equations 
will be the subject of future research. 

One of the principal motivations for developing the low Mach number equations and our numer- 
ical implementation was to model recent experiments on the development of giant concentration 
fiuctuations in the presence of sharp concentration gradients. We first studied giant fluctuations 
in a time-independent or static setting, as observed experimentally by inducing a constant concen- 
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tration gradient via a constant applied temperature gradient. Our simulations show that under 
conditions employed in experimental studies of the diffusive mixing of water and glycerol, it is 
reasonable to employ the Boussinesq approximation. At the same time, the results indicate that 
it is not accurate to employ a constant-transport-coefHcient approximation that is commonly used 
in theoretical calculations. 

We continued our study of giant concentration fluctuations by simulating the temporal evolution 
of a rough diffusive interface during the diffusive mixing of hard disk fluids. Comparison between 
computationally-intensive event-driven molecular dynamics simulations and our hydrodynamic cal- 
culations demonstrated that the low Mach number equations of fluctuating hydrodynamics provide 
an accurate coarse-grained model of fluid mixing. Special care must be exercised, however, in choos- 
ing the bare transport coefficients, especially the concentration diffusion coefficient, as these are 
renormalized by the fluctuations and can be strongly grid-dependent. Some questions remain about 
how to define and measure the bare transport coefficients from microscopic simulations, but we 
show that simply comparing particle and hydrodynamic calculations at large scales is a robust 
technique. 

The strong coupling between velocity fluctuations and diffusive transport means that determinis- 
tic models have limited utility at mesoscopic scales, and even macroscopic scales in two-dimensions. 
This implies that standard fluorescent techniques for measuring diffusion coefficients may not in 
fact be measuring material constants but rather geometry-dependent values. Fluctuating hydro- 
dynamic simulations of typical experimental simulations, however, are still out of reach due to the 
very large separation of time scales between mass and momentum diffusion. Surpassing this limi- 
tation requires the development of a semi-implicit temporal discretization that is stable for large 
time steps. Furthermore, it is also necessary to develop novel mathematical models and algorithms 
that are not only stable but also accurate in the presence of such large separation of scales. This 
is a nontrivial challenge if thermal fluctuations are to be included consistently, and will be the 
subject of future research. 
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Appendix 
Appendix A: Spatial Filtering 

In our spatial discretization, we chose to use centered differencing for the advective terms because 
this leads to a skew-adjoint discretization of advection [56J that maintains discrete fluctuation- 
dissipation balance in the spatially-discretized stochastic equations [4H . It is well-known that 
centered discretizations of advection do not preserve monotonicity properties of the underlying 
PDEs in the deterministic setting, unlike one-sided (upwinded) discretizations. Therefore, our 
spatio-temporal discretization can lead to unphysical oscillations of the concentration and density 
in cases where the cell Peclet number Pe = Ax ||i'|| /x is large. 

In the deterministic setting, Pe can always be decreased by reducing Ax and resolving the fine 
scale dissipative features of the flow. However, in the stochastic setting, the magnitude of the 
fluctuating velocities at equilibrium is 

where A^ is the volume of the hydrodynamic cell. Therefore, in two dimensions we have that the 
characteristic advection velocity magnitude is \\v\\ ^ Ax~^ This means that in two dimensions Pe 
is independent of the grid size and reducing Ax cannot fix problems that may arise due to a large 



cell Peclet number. For some of the simulations reported in Section VI we have found it necessary 
to implement a spatial filtering procedure that aims to reduce the magnitude of the fluctuating 
velocities while preserving their spectrum as well as possible at small wavenumbers. 

The ffitering procedure consists of applying a local averaging operation to the spatially- 
discretized random fields W and VV independently along each Cartesian direction. This local 
averaging smooths the random forcing and thus reduces the spectrum of the random forcing at 
larger wavenumbers. The specific filters we use are taken from Ref. [70]. For stencil width wp = 2, 
filtering a discrete field W in one dimension takes the form 

8 4 16 
In Fourier space, for discrete wavenumber Ak ~ kAx this local averaging multiplies the spectrum 
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oi W hy T (Afc) = 1 + O (Ak^) and therefore maintains the second-order accuracy of the spatial 
discretization. At the same time, the filtering reduces the variance of the fluctuating fields by 
about a factor of two in one dimension (a larger factor in two dimensions). The spectrum of the 
fluctuations can be preserved even more accurately if a stencil of width wj? = 4 is used for the local 

averaging, 

oQ 7 7 1 1 

giving a sixth-order accurate filter J^{Ak) = l + O {Ak^) and a reduction of the variance by about 
a third in one dimension. In two and three dimensions the filtering operators are simple tensor 
products of one-dimensional filtering operators. Note that we only use these filters with periodic 
boundary conditions. One can, of course, also use Fourier transform techniques to filter out high 

frequency components from the stochastic mass and momentum fluxes. 
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